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Abstract:  An  efficient  algorithm 
depth-separated  wwve  equation  in  g( 
media.  The  algorithm  has  been  bnilt 
pater  codes  called  SAFARI.  The  pa< 


es  been  developed  for  solving  the 
sral  fluid /solid  horisontally  stratified 
to  a  general- purpose  package  of  com- 
kge  consists  of  three  modules  provid¬ 


ing  plane  were  reflection  coefficients,  ftewHransmisaion  losses,  and  broadband 
pulse  response.  This  document  <$np  describes  the  mathematical  model  for 
seismo-acoustic  propagation  in  stratified  media.  Then  the  numerical  solu¬ 
tion  technique  is  outlined  followed  by  a  description  of  the  three  different 
SAFARI  modules  and  their  implementation.  ‘'Fiaally^t^e  actual  use  of  the 
different  modules  is  described,  including  a  detailed  discussion  on  the  numeri¬ 
cal  considerations  that  are  crucial  for  successful  use  of  this  type  of  numerical 
model.  SAFARI  is  applicable  to  a  wide  range  of  problems  in  many  disci¬ 
plines,  from  seismology  to  ultrasonics.  Here  its  use  is  illustrated  by  a  series 
of  examples  from  underwater  acoustics.  ^ 
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1.  Introduction 


The  seismo-acoustic  environments  encountered  in  nature  are  in  general  very  com¬ 
plex.  Sound  velocity  in  the  ocean  varies  with  both  range  and  depth  due  to  changes 
in  temperature,  salinity  and  density,  induced  by  solar  heating,  gravity,  currents  and 
eddies.  The  water  depth  can  vary  significantly,  in  particular  near  the  shore  and  at 
the  continental  rise.  In  addition,  the  bottom  will  often  be  characterised  by  irreg¬ 
ular  stratification  and  anisotropy.  Since  the  seismo-acoustic  propagation  is  highly 
dependent  on  the  material  properties,  an  exact  computer  modelling  would  require 
exact  knowledge  of  the  environmental  properties  to  the  smallest  detail,  which  is 
of  course  prohibitive.  Some  kind  of  approximation  therefore  has  to  be  applied  in 
order  to  obtsun  a  physical  model  of  the  environment  for  which  the  wave  equation 
can  be  solved  numerically.  Solution  techniques  such  as  the  finite-difference  tech¬ 
nique  [1]  and  the  finite-element  method  [2]  require  only  a  few  approximations  of  the 
environment.  These  techniques  are  therefore  the  most  general,  but  the  computa¬ 
tional  requirements  are  extensive,  in  reality  prohibiting  their  application  to  large 
scale  propagation  modelling.  None  of  the  more  efficient  techniques  are  applicable 
to  the  general  problem  since  they  are  based  on  r  *ific  assumptions  concerning  the 
environment.  The  parabolic  equation  technique  [4j  can  treat  propagation  in  a  range- 
dependent  environment,  but  the  shear  properties  of  the  ocean  bottom  are  ignored 
and  the  results  are  only  accurate  for  moderate  grazing  angles.  One  class  of  solution 
techniques  requires  the  environment  to  be  described  by  a  physical  model  for  which 
the  wave  equation  is  separable.  This  class  includes  the  normal  mode  [4]  technique, 
which,  however,  in  most  implementations  ignores  shear  and  is  limited  to  propagation 
at  grazing  angles  less  than  critical  (by  means  of  mode  coupling,  the  normal  mode 
technique  may  be  applied  to  non-seperable  problems  as  well  [5]).  Another  solution 
technique  requiring  the  wave  equation  to  be  separable  is  the  Fast  Field  (FFP)  or 
full  wavefield  technique  applied  in  SAFARI.  This  technique  yields  -  at  least  in  prin¬ 
ciple  -  an  exact  solution  to  the  wave  equation  in  a  horizontally  stratified  fluid/solid 
environment. 


The  principle  of  wave  equation  separation  for  horizontally  stratified  media  was  in¬ 
troduced  in  underwater  acoustics  by  Pekeris  [6j,  who  treated  t’j  problem  of  acoustic 
propagation  in  plane  layered  waveguides  using  simple  two-  and  three-layered  envi¬ 
ronmental  models.  Later,  Jardetsky  [7]  and  Ewing,  Jardetzky  and  Press  [8]  used  the 
same  technique  to  investigate  seismic  propagation  in  few-layered  waveguides.  The 
technique  was  to  apply  a  series  of  integral  transforms  to  the  Helmholtz  wave  equa¬ 
tion  to  reduce  the  origin^  four- dimensional  partial  differential  equation  (3  space 
dimensions  and  1  of  time)  to  a  series  of  ordinary  differential  equations  in  the  depth 
coordinate.  These  differential  equations  were  then  solved  analytically  within  each 
layer  in  terms  of  unknown  amplitudes  which  were  determined  by  matching  of  the 
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boundary  conditions  at  the  interfaces.  The  displacement  and  stress  quantities  were 
finally  determined  by  evaluation  of  the  inverse  integral  transforms. 

For  the  few-layered  cases  originally  presented,  the  linear  system  of  equations  in 
the  unknown  wavefield  amplitudes,  expressing  the  boundary  conditions,  can  easily 
be  solved  analytically.  For  more  complicated  environmental  models,  however,  this 
procedure  is  inconvenient  and  has  to  be  replaced  by  a  numerical  technique. 

Traditionally,  computation  of  the  depth  dependence  of  the  field  has  been  per¬ 
formed  by  means  of  propagator  matrix  methods  as  introduced  by  Thomson  [9]  and 
Haskell  [10].  The  propagator  matrix  approach  has  the  computational  advantage  that 
it  only  requires  a  small  amount  of  computer  memory  due  to  its  recursive  nature. 
However,  it  was  realised  quite  early  that  special  numerical  treatment  is  required  in 
order  to  ensure  numerical  stability,  and  several  modified  propagator  matrix  s diemes 
have  been  proposed.  In  general,  however,  these  have  resulted  in  much  more  time 
consuming  codes.  Further,  the  propagator  technique  is  not  well  suited  to  problems 
where  the  field  n as  to  be  determined  at  more  than  a  single  receiver  depth.  A  review 
of  the  propagator  approaches  is  given  by  Kenneth  [11]  who  himself  introduced  the  el¬ 
egant  invariant  embedding  formulation  [12]  which  has  the  interpretational  advantage 
that  arrivals  resulting  from  reflections  from  a  single  interface  can  be  isolated. 

The  propagator  matrix  approach  has  formed  the  basis  of  several  application  codes  in 
both  underwater  acoustics  [13]  and  seismology  [14].  The  so-called  Fast  Field  program 
developed  by  DiNapoli  [15]  applies  a  very  elegant  recursive  technique  to  determine 
the  depth- dependent  solution  for  many  horisontal  wavenumbers  simultaneously,  and 
is  therefore  extremely  efficient.  In  contrast  to  the  other  techniques,  however,  the 
depth-dependent  solution  is  approximate,  and  the  technique  is  only  applicable  to  a 
limited  class  of  fluid  problems. 

In  SAFARI  a  direct,  global  matrix  approach  is  taken  to  determine  the  depth  de¬ 
pendence  of  the  field  solution,  known  as  the  depth-dependent  Green’s  function. 
The  SAFARI  technique  is  in  fact  a  general  numerical  implementation  of  the  original 
solution  technique  of  Ewing  et  al.  [8],  but  implemented  using  efficient  numerical 
techniques  adopted  from  modern  finite-element  programs. 

In  what  might  be  termed  a  direct  global  matrix  or  ‘finite  wave  element’  approach, 
the  wavefield  in  each  layer  is  considered  as  a  superposition  of  the  field  produced  by 
an  arbitrary  number  of  sources  and  an  unknown  field  satisfying  homogeneous  wave 
equations.  These  unknown  fields  are  then  determined  from  the  boundary  condi¬ 
tions  to  be  satisfied  simultaneously  at  all  interfaces.  The  local  boundary  conditions 
at  each  interface  lead  to  a  linear  system  of  equations  in  the  Hankel  transforms  of 
the  potentials  in  the  adjacent  layers.  These  local  systems  of  equations  are  then 
straightforwardly  assembled  in  a  global  system  of  equations  expressing  the  bound¬ 
ary  conditions  at  all  interfaces  simultaneously.  The  resulting  global  coefficient  ma¬ 
trix  is  organised  to  be  block- bidiagonal  and  diagonally  dominant  in  close  analogy 
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to  the  global  stiffness  matrix  arising  in  the  finite-element  method.  Uncondition¬ 
ally  numerically-stable  solutions  are  determined  efficiently  by  gaussiaa  elimination, 
yielding  the  field  in  all  layers  simultaneously. 

Despite  the  analytical  equivtlence  of  the  propagator  and  global  matrix  solutions  for 
the  depth-dependent  Green’s  function,  there  are  a  number  of  important  advantages 
of  the  latter  technique  for  application  to  general  wave  propagation  problems  in 
underwater  acoustics  and  seismology  as  well  as  in  ultrasonics: 

e  The  field  produced  by  multiple  sources,  e.g.  in  phased  arrays,  is  easily  deter¬ 
mined  since  the  fields  produced  by  the  individual  sources  within  a  layer  are 
simply  superimposed,  and  no  dummy  interfaces  are  required  at  the  source 
depths. 

•  The  field  can  be  determined  at  any  number  of  receiver  depths  in  a  single  solu¬ 
tion  pass,  since  the  wavefield  potentials  are  found  in  all  layers  simultaneously. 

•  In  contrast  to  the  situation  for  techniques  based  on  propagator  matrices, 
mixed  problems  with  fluid,  solid  and  vacuum  layers  are  readily  treated  in  an 
efficient  manner,  as  the  global  system  of  equations  is  set  up  to  involve  only 
non- vanishing  potential  functions. 

•  Time  consuming  stability  assurance  problems  do  not  arise,  because  they  are 
removed  automatically  by  choosing  an  appropriate  coordinate  system  within 
each  layer  together  with  a  specific  organisation  of  the  global  system  of  equa¬ 
tions. 

•  The  solution  is  unconditionally  numerically-stable  for  any  number  of  layers. 

•  Due  to  the  global  nature,  many  operation?  can  be  vectorised,  making  the 
code  efficient  on  modern  array  and  vector  processors,  in  particular  in  cases 
with  many  layers. 

These  properties  make  the  global  matrix  approach  well  suited  for  implementation 
in  a  general  applications  program  like  SAFARI,  rendering  it  applicable  to  a  large 
range  of  seismo- acoustic  propagation  problems  in  underwater  acoustics,  seismology 
and  ultrasonics. 

Although  the  solution  for  the  depth-dependence  of  the  field  is  exact  to  within  ma¬ 
chine  accuracy,  it  is  a  common  characteristic  of  all  numerical  models  of  the  full 
wavefield  type  that  they  are  not  easy  to  use.  They  are  what  might  be  termed  'sci¬ 
entific'  numerical  models  which  require  the  user  to  possess  a  significant  knowledge 
of  the  physics  involved  in  wave  propagat;on  phenomena  in  stratified  waveguides  and 
at  the  same  time  be  confident  with  numerical  analysis.  This  is  due  to  the  fact  that 
the  evaluation  of  the  inverse  integral  transforms  has  to  be  made  numerically,  which 
requires  both  truncation  and  discretisation  of  the  integration  interval,  both  of  which 
introduce  errors  and  numerical  artifacts  which  have  to  be  reduced  to  insignificance. 
There  are  no  general  rules  to  follow,  because  the  necessary  sampling  depends  on  the 
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characteristics  of  the  actual  problem.  This  is  particularly  true  in  cases  where  only 
selected  parts  of  the  solution  are  of  interest.  The  numerical  integration  therefore 
cannot  be  fully  automated,  but  has  to  be  controlled  by  the  user. 

The  present  document  is  not  intended  for  teaching  waveguide  propagation  and  nu¬ 
merical  analysis  to  the  new  SAFARI  user.  Those  topics  are  treated  in  many  different 
textbooks.  Thus  the  intention  is  solely  to  present  a  self-contained  user's  guide  which 
should  introduce  the  new  user  to  the  SAFARI  model  and  at  the  same  time  act  as  a 
reference  manual  for  the  experienced  user. 

In  order  to  describe  the  basic  assumptions  underlying  the  SAFARI  codes,  the  math¬ 
ematical  model  for  wave  propagation  in  stratified  media  is  first  outlined.  Secondly, 
the  numerical  implementation  of  the  mathematical  model  is  discussed  in  detail,  pre¬ 
senting  first  the  global  matrix  solution  technique  particular  to  SAFARI  and  then  the 
different  numerical  integration  schemes  available  for  evaluating  the  total  field.  In 
the  third  part,  the  implementation  of  the  three  basic  SAFARI  modules  is  discussed 
in  general  terms;  the  precise  pr  edure  of  course  depends  on  the  actual  installation. 
Finally,  the  last  and  most  important  part  treats  the  process  of  running  the  codes, 
including  preparation  of  input  files  and  the  necessary  numerical  considerations  con¬ 
cerning  the  numerical  integrations  in  particular.  The  actual  use  is  illustrated  by 
a  series  of  characteristic  numerical  examples  from  underwater  acoustics.  Applica¬ 
tions  of  SAFARI  to  other  problems,  both  from  seismology  and  ultrasonics,  are  given 
in  [16-27]. 
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2.  The  environmental  model 

The  mathematical  description  of  the  seiamo- acoustic  propagation  used  in  SAFARI 
requires  the  environment  to  be  represented  as  a  horisontally  stratified  medium  as 
illustrated  in  Fig.  1.  Ail  interfaces  are  plane  and  parallel,  and  the  layer  properties 
are  range-independent.  The  layers,  including  the  upper  and  lower  halfspaces,  may 
be  either  viscoelastic  solids  or  fluids,  or  they  may  be  empty  space,  i.e.  vacuum. 


Interface  X 
Interface  £ 


Interface  m 


Layer  1:  Upper 

halfepacc 

Layer  2 

Layer  m 
Layer  m+1 


Interface  N-I 


layer  N:  Lower 

halfapace 


Fig.  J.  Horisontally  stratified  environment. 


The  solid  layers  are  required  to  be  homogeneous  and  isotropic  with  Lam6  constants 
X  and  fi  and  density  p.  The  corresponding  congressional  wave  speed  is 

Ce  =  n/(A  +  2  n)/p,  (1) 

and  the  shear  speed  is  _ 

e.  ^  \IWp-  (i) 


In  order  for  .the  solid  medium  to  have  positive  compressibility,  it  is  required  that  the 
bulk  modulus,  K  =  A  +  f/x,  be  positive.  A  physically  realistic  material  therefore 
requires  _ 

c,  <  \/0J5  cc.  (3) 
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In  the  fluid  layer*  the  density  p  must  be  constant,  whereas  the  nonvanishing  L&me 
constant  A  is  allowed  to  vary  with  depth  according  to  the  linear  law: 

t-t-t  =  az  +  6  (4) 

A(z) 

which  corresponds  to  the  following  depth  dependence  for  the  sound  speed.  Eq.  (1), 

ce(z)  =  s/l/(p{az  +  6)).  (5) 

The  special  caie  of  a  homogeneous  fluid  medium  corresponds  to  a  =  0. 
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3.  The  mathematical  model 

Only  two-dimensional  propagation  problems  can  be  treated  by  the  standard  SAFARI 
package,  although  a  special  three-dimensional  version  has  been  developed  [28].  Thus 
the  seismo- acoustic  field  has  to  be  either  plane  or  axisymmetric,  restricting  the 
source*  to  be  either  line  sources  perpendicular  to  the  plane  of  propagation  or  omni¬ 
directional  point  sources  placed  on  a  common  vertical  suds.  In  underwater  acoustics 
the  last  type  is  the  most  common,  and  therefore  we  here  describe  the  mathematical 
model  for  the  axisymmetric  case.  The  corresponding  model  for  the  plane  case  is 
described  in  [21]. 

The  full  wavefield  solution  technique  is  based  on  the  fact  that  for  a  horizontally 
stratified  environment,  it  is  possible  to  obtain  exact  integral  representations  for  the 
field  within  each  layer  in  terms  of  a  set  of  unknown  coefficients.  These  coefficients 
are  found  by  matching  the  boundary  conditions  simultaneuosly  at  all  interfaces,  and 
the  total  field  is  determined  by  evaluation  of  the  integral  representations.  In  this 
section  the  basic  principle  of  the  solution  of  the  wave  equation  by  depth-separation 
is  first  described.  Then  the  resulting  field  representations  will  be  given  for  both  fluid 
and  solid  media,  and  finally  the  pertinent  boundary  conditions  are  discussed.  The 
actual  numerical  implementation  is  discussed  in  Sect.  4. 


3.1.  THE  DEPTH-SEPARATED  WAVE  EQUATION 

A  cylindrical  coordinate  system  {r,  $,  z)  is  introduced,  Fig.  1,  with  the  z-axis  pass¬ 
ing  through  the  sources,  making  the  field  independent  of  azimuthal  angle  0.  For 
the  isotropic  media  considered  here,  the  seismo- acoustic  wavefield  can  then  be  ex¬ 
pressed  in  terms  of  scalar  wavefield  potentials  ¥(r,  z,  t),  which  satisfy  the  linear  wave 
equation 

<•> 

Here  c(r,  x)  is  a  wave  speed  and  F,(r,  z,  t)  a  forcing  term  which  in  the  present  case 
is  due  to  the  seismo-acoustic  sources.  We  now  utilize  the  Fourier  transform  pair 


f(u>)  = 

(7) 

fCO 

F(t)  = 

/  /(u/)e,wt  du>. 

J — oo 

(8) 
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By  applying  the  forward  transform,  Eq.  (7),  to  the  wave  equation,  Eq.  (6),  the 
frequency  domain  wave  equation  is  obtained: 

(V*  4-  *£,(r,  *))  *(»•.  w)  =  /*(r,  z,u>),  (9) 

where  u>  is  the  angular  frequency  and  km(r,  z)  is  the  medium  wavenumber, 

km(r,z)  =  u>/c(r,z).  (10) 

Since  the  environment  is  range  independent  and  the  sources  are  confined  to  the 
vertical  axis,  Eq.  (9)  can  be  rewritten  as 


(V*  +  tl(z))  (11) 

In  the  following  the  u  dependence  will  be  generally  suppressed. 

Next,  we  employ  the  Hankel  transform  pair 


9(k)  = 

f  G(r)J0(kr)rdr , 

Jo 

(12) 

G(r)  = 

f  g(k)J0{kr)kdk, 

Jo 

(13) 

where  k  is  the  horizontal  wavenumber,  and  apply  the  forward  transform,  Eq.  (12), 
to  Eq.  (11),  obtaining  the  depth-separated  wave  equation 

(■|r> '  *<*■ 2) = (14) 

i.e.  an  ordinary  differential  equation  in  depth.  The  solution  is  the  sum  of  a  particular 
solution  9(k,z)  to  Eq.  (14)  and  any  linear  combination  of  the  two  independent 
solutions  9~(k,  z)  and  9+(k,  z)  to  the  homogeneous  equation 

(•^r *(*.*)  =  °-  (is) 

The  total  i  lution  for  the  depth  dependence  of  the  field,  the  so-called  depth-depen¬ 
dent  Green’s  function,  is  therefore 

»(*,  *)  =  *(*,  *)  +  i4-(ft)«"(fc,  z)  +  A+(k)9+(k,  z),  (16) 

where  i4~(ft)  and  j4+(k)  cue  arbitrary  coefficients  to  be  determined  from  the  bound¬ 
ary  conditions.  The  particular  solution  to  Eq.  (14)  is  most  conveniently  chosen  to 
be  the  field  produced  by  the  sources  in  the  absence  of  boundaries. 
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When  the  unknown  coefficients  are  found,  the  total  field  at  the  angular  frequency  u> 
is  found  at  any  range  r  by  carrying  out  the  inverse  Hankel  transform,  Eq.  (13).  Sim¬ 
ilarly,  the  time  response  is  determined  by  evaluating  the  inverse  Fourier  transform, 
Eq.  (8). 

The  basic  property  of  the  full  wavefield  solution  technique  is  to  restrict  the  depth 
dependence  of  t*,(x)  to  cases  where  Eqs.  (14)  and  (15)  can  be  solved  analytically, 
limiting  the  numerical  effort  to  determining  the  unknown  coefficients  A~(k)  and 
A+(k)  from  the  boundary  conditions  and  to  evaluating  the  inverse  integral  trans¬ 
forms. 

In  the  following,  the  analytical  field  representations  will  be  given  for  the  media 
included  in  the  SAFARI  model,  followed  by  e,  discussion  of  the  boundary  conditions 
to  be  satisfied  at  the  interfaces. 


8.J.  HOMOGENEOUS  FLUID  MEDIUM 


Fbr  an  ideal  fluid,  the  equation  of  motion  for  angular  frequency  <*>  is  easily  shown 
to  be  satisfied  if  the  displacements  are  expressed  in  terms  of  a  scalar  displacement 
potential  as 

d* 

U~  dr 
6$ 
w~  dz 

where  u  and  w  are  the  radial  and  vertical  components,  respectively,  and  the  potential 
#(r,  z)  satisfies  the  wave  equation,  Eq.  (9).  The  depth  dependence  of  the  field  is 
therefore  of  the  form  given  by  Eq.  (16). 

In  the  case  of  a  homogeneous  fluid  layer,  the  medium  wavenumber  is 


(17) 

(18) 


^m(*)  —  hm»  (19) 

i.e.  the  compressions!  velocity  ce  is  constant  not  only  in  range  but  also  in  depth,  and 
the  homogeneous  solutions  obtained  from  Eq.  (15)  are  simply  exponential  functions, 

i+(k,z)  =  eat,  (20) 

{k,z)=  e-**,  (21) 

where 

<*'*)  =  -  hm-  (22) 
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If  no  sources  are  present  in  the  layer,  the  total  field  then  follows  directly  from  Eq.  ( 13) 

[. A~e~a '  +  A+eat]  J0{kr)kdk.  (23) 

In  physical  terms,  Eq.  (23)  is  a  decomposition  of  the  total  wavefield  in  upgoing  (e°*) 
and  downgoing  (e-a*)  conical  waves  with  horizontal  vavenumber  k. 

The  boundary  conditions  for  fluid  layers  ,  discussed  in  Sect.  3.5,  involve  the  vertical 
displacement  w  and  the  normal  stress  <rxt  (equivalent  to  negative  pressure  in  a  fluid 
medium).  The  vertical  displacement  follows  from  Eq.  (18), 

«jfr,2)=  r  \-aA~e-az  +  aA+eat]Mkr)kdk,  (24) 

Jo 

whereas  the  normal  stress  follow  from  Hooke’s  law, 

c„(r, z)  =  AVJ$(r, *1  =  -pcv*$(r, z)  =  -pJ1  f  \A~c~a*  +  A+ea~]  J0(kr)kdk. 

Jo 

(25) 


If  a  source  is  present  in  the  layer,  the  particular  solution  to  Eq.  (14)  has  to  added. 
In  the  case  of  an  omnidirectional  point  source,  the  term  /,(z,w)  in  Eq.  (11)  takes 
the  form 

/.(z.w)  -  -Sw$(r  -  z, ),  (26) 

where  Su  is  the  source  strength  and  z,  is  the  source  depth.  It  can  be  shown  [30], 
that  in  this  case  the  particular  solution  to  Eq.  (14)  is 


•(*.*) 


Swe-*1*—1 
4ir  a 


(27) 


and  the  corresponding  source  contribution  to  the  total  field  again  follows  from 
Eq.  (13), 

*(r,z)=- fj0  - a - Jo^r)kdk.  (28) 

If  more  sources  are  present  in  the  layer  their  contributions  are  simply  superimposed. 
The  displacements  and  stresses  are  derived  analogously  to  Eqs.  (24)  and  (25)  and 
become 


tD(r,z)=-7!:i  /  sign(z  -  z,)e  Q|*  **l/o(kr)kdlf,  (29) 

4*  J o 

S.,ow 1  e~a I*-1*! 

<7**(»\2)= - -  /  - J0(kr)kdk.  (30) 

4  t  y0  a 


-10- 


Saolantobn  sr-hs 


3  -  The  mathematical  model 


U.  INHOMOGENEOUS  FLUID  MEDIUM 

la  a  real  ocean  environment  the  sound  velocity  varies  with  depth,  and  it  is  well 
known  that  this  variation  as  significant  impact  on  the  acoustic  propagation.  The 
ocean  waveguide  can  therefore  not  in  general  be  represented  by  a  single  homogeneous 
fluid  layer.  It  can  be  shown,  however,  that  if  the  range-independent  ocean  waveguide 
is  represented  by  an  increasing  number  of  homogeneous  layers,  a  numerical  solution 
based  on  the  Add  representation  for  homogeneous  layers  will  converge  towards  the 
comet  solution.  However,  a  satisfactory  convergence  requires  the  layers  to  be  less 
than  one  quartet  of  a  wavelength  thick.  Such  a  technique  is  therefore  only  convenient 
for  low  frequency  propagation  in  shallow  water,  and  modelling  of  deep  water  or  high 
frequency  propagation  does  require  some  sort  of  velocity  profile  to  be  incorporated. 

Since  it  is  still  necessary  that  Eq.  (14)  can  be  solved  analytically  it  is  obvious  that  the 
ocean  cannot  in  general  be  represented  by  a  single  layer,  but  has  to  be  approximated 
by  a  series  of  layers,  within  which  the  depth  dependence  of  the  field  has  an  analytic 
representation.  A  few  examples  of  sound  speed  interpolation  functions  for  which  this 
is  possible  are  given  in  [15].  However,  the  actual  choice  of  interpolation  function  is 
not  very  critical  since  the  profile  is  usually  measured  at  discrete  depths  with  a  finite 
uncertainty.  In  the  SAFARI  code  the  solution  corresponding  to  the  sound-speed 
variation  given  in  Eq.  (5)  has  been  implemented.  In  this  case  the  homogeneous 
depth-separated  wave  equation,  Eq.  (15),  becomes 

(■J?  "  + l»)  *(*■  *)  =  o.  (31) 


By  introducing  the  variable  transformation, 

c  =  c  '*/*(«  +  d) 

=  (. Va)_,/,(*J  -  pu >*{az  +  h)),  (32) 

the  following  equation  is  obtained: 

(•|r  -  <)  *<0  =  0.  (33) 

This  is  a  special  form  of  the  Bessel  differential  equation,  for  which  two  independent 
solutions  are  the  Airy  functions  Ai(()  and  £t((),  [29].  Independent  homogeneous 
solutions  to  Eq.  (31)  are  therefore 

i+(ktz)  =  Ai(c’''*(k*-(cz  +  d))t 
*-(*,  z)  =  Bi( c-*/*(fc*  -  («  +  d)). 
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The  total  field  at  angular  frequency  u  is  again  obtained  by  evaluation  of  the  inverse 
Hankel  transform,  Eq.  (13), 


*(r,r)  =  [“[A+AiiO  +  A~  Di(())Mkr)kdk, 

Jo 

with  the  corresponding  displacement  and  stress  given  by 

w(r,z)  =  -c1'*  JJ°[A+Ai,(C)  +  A~  Bi'iOW^k&k, 

<r„(r:r)=  -po,1  ["[A*  Ai{Q  +  A~  Bi(Q]Mkr)k&k, 
Jo 

where  a  prime  denotes  the  derivative  with  respect  to  the  argument. 


(36) 

(37) 

(38) 


When  a  source  is  present  in  the  layer  the  corresponding  particular  solution  to 
Eq.  (14)  again  has  to  be  added.  However,  even  in  the  case  of  a  point  source  it  is 
not  straightforward  to  find  such  a  particular  solution.  We  will  apply  a  small  ’trick’ 
to  obtain  analytical  representations  for  the  depth  dependence  of  the  field  produced 
by  a  point  source  at  depth  z,  in  an  infinite  medium  with  the  sound  velocity  given 
by  Eq.  (5).  A  thin  homogeneous  layer  of  thickness  e  with  sound  speed  cc  =  cc(rt) 
is  introduced  containing  the  source,  Fig.  2.  In  the  limit  e  ->  0,  the  solution  of  this 
layered  problem  converges  to  the  solution  of  the  original  problem. 


Fig.  2.  Point  source  in  inhomogeneous  fluid  medium. 
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This  S  -Ur«/  problem  is  solved  in  exactly  the  same  way  in  which  the  general  multi¬ 
layered  problem  is  solved  in  SAFARI.  As  an  illustration  of  this  solution  technique, 
we  shall  here  go  through  the  details.  In  Fig.  2  the  sound  speed  variation  is  shown 
few  the  case  in  which  a  >  0,  i.e.  lim,_oo  ce(z)  =  0  and  limA__s/a  ce(z)  =  oo.  The 
argument  (  to  the  Airy  functions,  Eq.  (32),  therefore  has  the  limits 


lim  C  =  -oo. 

(39) 

•  -♦00 

lim  C  =  oo- 

— s-OO 

(40) 

Since  no  sources  are  present  in  the  two  halfspaces,  the  field  in  both  is  given  by 
Eq.  (36).  The  field,  however,  must  satisfy  the  radiation  condition  for  z  — *  ±oo, 
which  immediately  reduces  the  number  of  unknown  coefficients.  Since  the  field 
must  be  limited  and  lim<_0o  Bs(C)  —  oo,  the  depth  dependence  of  the  field  in  the 
upper  halfspace  must  be  of  the  form 

*,(M)  =  ^At«).  (41) 

Similarly,  the  radiation  condition  requires  that  only  downgoing  waves  exist  for  z  -+ 
oo.  Because  of  the  asymptotic  behaviour  of  the  Airy  functions  for  (  — ►  -oo  [29],  it 
is  required  that  the  field  in  the  lower  halfspace  be 

*,(k,z)  =  A,-(A»(C)  +  <Bs(0).  (42) 


For  the  intermediate  isovelocity  layer,  the  depth  dependence  of  the  field  is  directly 
obtained  from  Eqs.  (20),  (21)  and  (27)  as 


$,(*,  z)  =  Afe~aM  +  Afeat  + 


Sv 

4r  a 


(43) 


The  next  step  is  to  satisfy  the  boundary  conditions  of  continuous  vertical  displace¬ 
ment  and  stress  at  the  two  interfaces,  which  are  so  close  together  that  they  can  both 
be  considered  to  be  at  depth  z,,  leading  to  the  following  system  of  equations: 

Atf-c'/’Ai'U.))  -  V(-A;  +  A})  =  |s,  (44) 

-A+Ai(C.)  +  (A;  VAi )  =  (45) 

o(-^r  +  A})-  A;(-c‘'’(Ai'( (.)  +  ifli'K,)))  =  (48) 

-MT  +  >if)  +  A;(Ai((.)  +  (47) 
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All  terms  involving  the  coefficient*  A,  and  A  f  ere  eeeily  eliminated  by  pairwise 
addition,  end  the  resulting  two  equations  give  the  following  solutions: 

.+  _  _£. _ 2c-'^(Ai((.)  +  iSi((.)> _ 

4<r  +  iSi((.))  -  -U((.)Mi'(C.)  +  ifli'K.))’  '  ’ 

_  _£» _ 2c-‘/Mi((.) _ 

*  _  4r  Ai-((.)(Ai((.)  +  +  iB.'«.))  '  *’ 


The  source  field  representations  for  the  case  where  the  sound  speed  increases  with 
depth,  i.e.  a  <  0,  can  of  course  be  determined  Li  the  same  way,  but  it  can  also  be 
determined  right  away  by  symmetry  considerations.  The  choice  of  depth  axis  x  is 
arbitrary,  and  we  can  therefore  perform  the  variable  transformation  z  —>  -z.  This 
will  also  change  the  sign  of  a  in  Eq.  (5),  and,  as  can  be  observed  from  Eq.  (32),  (  is 
then  invariant  to  this  transformation,  and  the  results  above  are  therefore  still  valid; 
they  just  have  to  be  interchanged  between  the  two  halfspaces.  Th''  held  produced 
by  a  source  in  an  inhomogeneous  fluid  layer  therefore  has  the  integral  representation 


$(r,s)  =  -^jf“j0(*r)td* 

_ 2c-*/»(A«K.)  +  »g»«,))Ai(C) _ 

Ai'(C.)Mi(C.)  +  iBi( C.))  -  A*(C,)(j4*'(C»)  +  »2K'(C.))’ 

<*(*  -  *.)  <  0, 

_ 2c~1^8Ai(C«)(At(C)  4-  iBi(Q) _ 

At'K.)(A»«.)  +  »!?*(£♦))  -  Ai(C«)(A*'(C«)  + 

<*(*  -  2.)  >  0. 


(50) 


The  resulting  stresses  and  displacements  are  derived  as  they  are  for  the  homogeneous 
case.  Note  that,  as  was  the  casejor  a  source  i  a  homogeneous  layer,  the  kernel 
in  the  integral  representation  for  $(r,  z)  is  continuous  at  the  source  depth,  but  the 
depth  derivative  is  discontinuous.  If  more  sources  are  present,  their  contributions 
are  simply  superimposed. 
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1.4.  HOMOGENEOUS  SOLID  MEDIUM 

Fbr  a  homogeneous  and  isotropic  elastic  solid,  it  can  be  shown  [8],  that  the  equation 
of  motion  is  satisfied  if  the  displacements  are  expressed  in  terms  of  scalar  potentials 

{*,*>« 


u(r,  z) 
u?(r,  z) 


(61) 

.  i  a  a ..  . 

(62) 

where  the  potential!  in  the  absence  of  sources  satisfy  the  uncoupled  homogeneous 
wave  equations 


(63) 

*("•*.  0  =  °.  (M) 

cc  and  c ,  being  the  compresaional  and  shear  velocities,  respectively,  given  in  terms  of 
the  Lam4  constants  in  Eqs.  (1)  and  (2).  Both  these  equations  are  of  the  same  form 
as  Eq.  (6)  and  can  therefore  be  depth-separated.  As  in  the  case  of  a  homogeneous 
fluid,  the  constant  wave  speeds  lead  to  simple  exponential  functions  in  depth,  and  the 
potentials  therefore  have  the  following  integral  representations  for  angular  frequency 
u>: 


#(r,  x)  =  j fV-e—*  +  A+ca*]/0(*r)fc  d  k,  (65) 

*(r,  z)  =  [“[B-e-**  +  B+e*9)Mkr)  d*.  (66) 

Jo 

where 

«(*)  =  (6)) 

m  =  (68) 

The  medium  wavenumbers  and  km  are  for  compressional  and  shear  waves,  re¬ 
spectively.  By  means  of  Eqs.  (51)  and  ,52),  the  following  integral  representations 
are  obtAned  for  the  dispiarerwuts, 

«(r,r)=  jT  [-iA-e-0*  -  bA+e°*  +  Jt{kr)kAk%  (59) 

*.  (r,x)  =  j r  [-aA'e"<"  +  aA+eaj  +  kB~e~fi9  +  kB+e*9)  J9{kr)kdkt  (60) 
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ana  the  normal  jtress  oMt  and  tangential  stress  <rTt  follow  from  Hooke's  law, 

,  v  _  .  dw  .  du 

*««(»•,  z)  =  (A  +  2/i)  ^ 

-2 kp[B-t-fi‘  -  S^1]]  >M*r  )Jfcdk, 

i  \  ( du  &w\ 

<r^z)=ti{^  +  ^) 

=  n  J~\2kct[A-e-a*  -  j4+ea*| 

-  (2**  -  kDiB-e'0'  +  B+  e^*]j  Jx(kr)k dk. 


In  a  solid  medium,  integral  representations  for  several  types  of  sources  can  be  de¬ 
rived,  including  point  forces  and  force  couples.  A  number  of  examples  are  given  by 
Ha  krider  [31].  In  the  SAFARI  code,  two  source  types  are  available  for  solid  media. 
One  is  a  congressional  point  source  equivalent  to  the  one  used  in  fluid  media.  The 
compressional  point  source  involves  the  compressional  potential  $  only,  with  the 


integral  representation 

*(r,*)  =  7^  /  - J0(kr)kdk.  (63) 

iw  J0  a 

The  corresponding  displacement  and  stresses  again  follow  from  Eqs.  (51), (52)  and 
Hooke's  law: 

S  ke  -<*!*—*•  I 

S(r,z)  =  - Jx(kr)kdk,  (64) 

4ir  y0  a 

tD(r,z)  =  jHsign^  -  zt)e-a'‘—'j0(kr)kdk,  (65) 

?..(■■.<)  =  («’  -  o — - — w>‘«,  (66> 

9„(r,x)=  ^  /“jt.ignfr-  • •lj,(*r)tdt.  (67) 

The  second  implemented  source  type  is  a  vertical  point  force,  with  the  following 
integral  representations  for  the  potentials  [8]: 

*(r,z)  =  J"  z.)e-a'‘~xAMkr)kdk,  (68) 
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t(r,  *)  =  £  J  t——j0(kr)kdk.  (89) 

Integral  representations  for  the  displacements  and  stresses  are  again  obtained  from 
Eqs.  (51), (52)  and  Hooke’s  law: 

«(r,,)=  jTsign(z  -  zt)k  [e-l*— '  -  J<(kr)kdky  (70) 

u?(r,  z)  -■  J”  [ae-0!*-1*1  -  k'fi-'e’* I*"**']  Jo(kr)kdk,  (71) 

5„(r,  z)  =  ^  J  sign(r  -  zt) 

x  [(2fc*  -  ~  Jo(kr)k dk,  (72) 

9r.(r,x)  =  ^  ^kae-**1*"**1  -  (*S/Tl  +  k/3)e~^'^]  Jx(kr)kdk.  (73) 


3.5.  BOUNDARY  CONDITIONS 

The  field  at  each  interface  now  has  two  distinct  integral  representations,  one  from 
the  layer  above  and  one  from  the  layer  below.  Depending  on  the  type  of  interface, 
a  certain  set  of  boundary  conditions  has  to  be  satisfied. 

•  If  thu  interface  is  separating  two  fluid  layers,  the  vertical  displacement  w  and 
the  normal  stress  <rM»  have  to  be  continuous.  If  one  of  the  media  is  a  vacuum, 
the  normal  stress  <r„  must  vanish. 

•  If  one  of  the  levers  is  solid  and  the  other  is  a  vacuum,  both  <ra«  and  <rri 
must  vanish,  whereas  the  fluid/solid  interface  requires  that  w  and  <r„  must 
be  continuous  and  <rTt  must  vanish. 

e  At  a  welded  interface  between  two  solid  media,  tu,  u,  <rMM  and  <rr*  all  have  to 
be  continuous. 

The  boundary  conditions  for  the  different  types  of  interfaces  are  summarised  in 
Table  1. 

Since  the  boundary  conditions  have  to  be  satisfied  at  all  ranges  r,  it  is  obvious  that 
they  must  be  satisfied  by  the  kernels  in  the  integral  representations  as  well.  The 
radiation  conditions  for  z  — *  ±oo  together  with  the  conditions  to  be  satisfied  at  all 
interfaces  simultaneously,  lead  to  a  linear  system  of  equations  in  the  unknown  kernel 
coefficients  A~  ,A+,B~  and  f?+.  In  principle,  this  system  has  to  be  solved  for  all 
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Table  1 

Boundary  conditions* 


Type 

Field  parameter 

w  u 

Or. 

(laid/ vacuum 

. 

0 

. 

fluid/fluid 

= 

- 

fluid/solid 

= 

= 

0 

solid/ vacuum 

- 

0 

0 

solid  /  solid 

= 

- 

*  Symbol*  used:  =,  continuous;  0,  vanishing;  not  in¬ 
volved. 


values  of  the  horizontal  wavenumber  Jfc,  and  the  total  field  can  then  be  determined 
by  evaluating  the  inverse  transforms.  Except  for  a  few  trivial  cases,  however,  both 
the  solution  of  the  linear  system  of  equations  and  the  evaluation  of  the  inverse 
transforms  have  to  be  done  numerically,  requiring  truncation  and  discretisation  of 
the  horisontal  wavenumber  axis. 


1.1  ATTENUATION 

It  is  well  known  that  sound  propagating  in  the  ocean  bottom  undergoes  a  significant 
attenuation  due  to  the  dissipation  of  seismo- acoustic  energy  into  heat.  It  is  there¬ 
fore  crucial  to  a  realistic  prediction  of  the  propagation  characteristics  that  volume 
attenuation  be  taken  into  account.  In  the  present  full  wavefield  solution  technique 
this  is  straightforwardly  done  as  shown  below. 

Assume  a  plane  harmonic  wave  of  angular  frequency  u>  propagating  in  a  homogeneous 
medium  along  the  positive  e-axis  of  a  cartesian  coordinate  system.  With  the  present 
choice  of  time-frequency  transform  pair,  Eqs.  (7), (8),  such  a  wave  has  the  form 

*•(*,0  =  >,  (74) 

where  km  is  the  medium  wavenumber  for  either  compressional  or  shear  waves  and 
A  is  the  amplitude.  If  im  is  real,  this  wave  has  constant  amplitude  for  sdl  ranges  c. 
It  can  be  shown  [32),  however,  that  viscoelastic  attenuation  csm  be  accounted  for  by 
letting  the  medium  wavenumber  km  be  complex,  i.e. 

km  —  km(l  d)i  >  0.  (76) 
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X  .  me  wavefteld  becomes 

F(*,  t)  =  Ae***"*  V*"*-*— >.  (76) 


Now  the  amplitude  is  decaying  exponentially  in  range  as  is  required  for  linear  vis¬ 
coelastic  fluid  and  solid  media.  Since  the  full  wavefield  solution  technique  is  based 
on  plane- wave  decomposition  as  described  above,  it  is  obvious  that  a  viscoelastic 
attenuation  is  taken  into  account  by  letting  the  medium  wavenumbers  and  therefore 
also  the  Lam4  constants  be  complex: 


A  =  A  +  iA\  (77) 

fi  =  f*  +  «>'  •  (78) 


It  has  been  experimentally  observed  [8,32],  that  most  solid  media  exhibit  an  atten¬ 
uation  that  increases  linearly  with  frequency,  i.e.  6  =  const.  For  these  solids, 


A*  +  2fi'  _  _1_ 
A  +  2/i  Qe' 

_ L. 

n  Q»' 


(79) 

(80) 


where  Qe  and  Q ,  are  constants.  By  inserting  the  complex  Lam6  constants  in  Eqa.  (1) 
and  (2),  we  get  from  Eqs.  (10)  and  (75)  the  following  value  of  6  for  compressions! 
and  shear  waves  respectively,  assuming  QC,Q,  >  1: 


6.= 


1 

2  Qe' 

(81) 

1 

2  q: 

(82) 

In  underwater  acoustics  it  is  more  common  to  express  the  linear  frequency- dependent 
attenuation  in  dB/A  (where  A  is  the  wavelength): 


7  =  -20  log 


g(*+jM) 


-20  log  [e  =  40*£loge  ~ 


(83) 


Although  the  solution  technique  is  valid  for  any  value  of  the  attenuations,  these 
must  be  physically  meaningful.  Thus  it  is  required  that  a  pure  dilatation  of  a  solid 
medium  does  not  produce  energy.  Therefore  the  bulk  modulus,  K  =  A  -f  must 
have  a  positive  imaginary  part,  which  is  easily  shown  to  require 


7. 

7e 


(84) 
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The  linear-elastic  fluid  media  are  merely  limiting  cases  of  solids  with  ft  -»  0,  and 
the  attenuation  is  therefore  introduced  accordingly. 

In  addition  to  the  volume  attenuation  described  above,  scattering  at  rough  inter¬ 
faces  also  attenuates  the  stress  waves  in  stratified  environments.  Since  the  depth- 
separation  solution  technique  requires  all  interfaces  to  be  plane  and  parallel,  a  deter¬ 
ministic  treatment  of  rough  surface  scattering  is  not  possible  by  means  of  the  present 
method.  Such  a  deterministic  treatment  is,  however,  only  of  academic  interest,  since 
the  roughness  appearing  in  nature  is  cither  rapidly  changing  (as  in  the  case  of  ocean 
surface  waves)  or  is  not  sufficiently  well  described  to  permit  direct  modelling.  The 
statistical  roughness  characteristics,  on  the  contrary,  are  relatively  slowly  varying 
and  can  be  more  easily  determined  even  for  the  ocean  bottom. 

If  the  roughness  is  small  compared  to  the  acoustic  wavelength  involved  and  the 
roughness  characteristics  are  stationary  in  both  space  and  time,  it  is  possible  by 
means  of  a  perturbational  approach  to  modify  the  boundary  conditions  in  such  a 
way  that  the  effect  of  the  roughness  on  the  coherent  (mean)  part  of  the  field  can  be 
taken  into  account  by  the  present  solution  technique.  Such  a  technique,  originally 
developed  for  fluid  media  only,  has  been  generalised  to  interfaces  between  solid 
media  and  implemented  in  the  SAFARI  code  [25].  The  actual  perturbation  theory  is 
not  of  importance  to  the  general  use  of  the  code,  and  therefore  only  the  underlying 
assumptions  concerning  the  roughness  characteristics  will  be  given  here.  For  details 
reference  is  made  to  [25]. 

The  interface  elevation  at  a  point  with  the  horizontal  coordinates  given  by  the  vector 
i'  is  y(r).  This  function  is  only  determined  in  terms  of  its  correlation  function  JV(R), 
defined  by 

=  <7(ri)7(rj))>  R  =  rx  -  r,  (86) 

or  in  terms  of  the  two-dimensional  Fourier  transform  of  the  interface-rou  hness 
power  spectrum  P(a), 

N(R)  =  <7*>^  J  ^aP(a)  e'u  R,  (86) 

where  (7*)  is  the  .ms  value  of  the  roughness.  In  SAFARI  only  isotropic  roughness 
can  be  treated,  i.e.  the  interface  correlation  depends  only  on  the  horizontal  distance 
|R|  between  two  points,  not  the  direction.  Further,  the  spectrum  -  of  course  also 
isotropic  -  is  assumed  to  be  gaussian: 


P(a)  =  L'e-^'W7,  (87) 

where  L  is  a  characteristic  correlation  length  for  the  interface  roughness.  The  inter¬ 
face  roughness  statistics  are  therefore  described  entirely  by  the  rms  roughness  (7*) 
and  the  correlation  length  L.  These  parameters  may,  of  course,  vary  from  interface 
to  interface.  For  the  perturbation  expansion  to  be  valid,  it  is  required  that  7/L  <  1. 


-20- 


Saclantcen  sr-;is 


4  -  Numerical  solution  technique 


4.  Numerical  solution  technique 


The  numerical  solution  of  the  full  wavefield  problem  divides  naturally  into  three 
parts.  First  the  depth-dependent  Green’s  function  is  found  at  a  discrete  number  of 
horizontal  wavenumbers  for  the  selected  receiver  depths.  Secondly,  the  wavenumber 
integral,  Eq.  (13),  is  evaluated,  yielding  the  transfer  function  at  the  selected  depths 
and  ranges.  Finally,  upon  repetition  of  the  first  two  steps  at  selected  frequencies, 
the  frequency  integration,  Eq.  (8),  is  performed  to  yield  the  total  response  in  time. 

It  is  clear  that  the  three  steps  are  not  independent.  In  particular,  the  frequency  and 
wavenumber  sampling  cannot  be  chosen  arbitrarily,  but  have  to  satisfy  the  require¬ 
ments  of  the  sampling  theorem.  The  frequency  interval  depends  on  the  bandwidth 
of  the  source  and  the  time  window  must  be  sufficiently  large  to  contain  the  whole 
signal  in  order  to  avoid  wrap-around,  which  again  imposes  requirements  on  the  fre¬ 
quency  sampling.  Similarly,  the  wavenumber  sampling  must  be  sufficiently  large  to 
allow  accurate  integration  of  the  often  rapidly  varying  and  oscillating  integrand  in 
Eq.  (13). 

Depending  on  the  time/frequency  and  range/depth  requirements,  the  determination 
of  the  depth-dependent  Green’s  function  may  have  to  be  performed  a  substantial 
number  of  times.  The  efficiency  of  the  code  is  therefore  highly  dependent  on  this 
part,  and  it  is  precisely  here  -  in  the  global  matrix  approach  -  that  the  SAFARI  code 
differs  in  computational  approach  from  recursive  propagator  matrix  techniques  used 
in  earlier  codes  of  the  same  type. 

The  global  matrix  solution  technique  forms  the  basis  for  all  three  SAFARI  modules. 
SAFARI-FIPR  calculates  plane-wave  reflection  coefficients  for  an  arbitrarily  strati¬ 
fied  halfspace.  SAFARI-FIP  first  calculates  the  depth-dependent  Green’s  function 
and  then  evaluates  the  wavenumber  integrals  to  obtain  the  single-frequency  seismo- 
acoustic  field  at  selected  receiver  depths  and  ranges.  The  last  module,  SAFARI-FIPP, 
performs  the  same  operation  for  a  number  of  frequencies,  and  after  weighting  with 
a  source  spectrum  it  determines  the  pulse  response  at  any  receiver  by  evaluating 
the  frequency  integrals.  In  this  section  we  will  describe  in  detail  the  numerical 
techniques  applied  in  all  three  program  modules. 
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4.1.  DEPTH-DEPENDENT  GREEN’S  FUNCTION 

The  use  of  digital  computers  for  numerically  solving  field  problems  in  continuous 
media  in  general  requires  some  kind  of  discretisation.  One  possibility  is  to  set  up 
the  exact  field  equations  for  the  continuum  and  subsequently  find  an  approximate 
solution  by  discretisation  of  the  equations.  This  is  the  approach  taken  in  finite 
difference  techniques.  The  other  possibility  is  to  discretise  the  medium  itself  and 
then  in  effect  let  the  computer  numerically  determine  an  exact  solution  to  this  now 
approximate  problem.  The  best-known  example  of  this  latter  approach  is  the  fi¬ 
nite  element  method  [2].  This  technique  is  based  on  the  division  of  the  continuum 
into  finite  blocks  or  elements  connected  to  each  other  at  a  finite  number  of  discrete 
nodes.  Exact  local  solutions  for  the  single  elements,  together  with  the  continuity 
between  the  elements,  concentrated  in  the  nodes,  lead  directly  to  an  exact  solu¬ 
tion  for  the  approximated  global  problem.  This,  however,  requires  the  solution  of 
very  large  linear  systems  of  equations  in  the  unknown  degrees  of  freedom  (typically 
node  displacements).  The  unavailability  of  sufficient  computing  power  until  the  last 
decade  delayed  the  widespread  use  of  this  powerful  technique,  which  today  is  the 
most  widely  employed  methodology  in  structural  and  fluid  mechanics,  for  example. 

In  this  regard,  it  is  easily  observed  [19,23]  that  the  integral  transform  solution  of 
the  wave  equation  for  horizontally  stratified  environments  is  of  the  ‘finite  element’ 
category.  The  local  expressions  for  the  depth-dependent  Green’s  function  are  exact 
within  each  layer  (element),  and  thus  an  exact  global  solution  can  be  obtained 
directly  from  the  boundary  conditions  to  be  satisfied  at  all  interfaces  (nodes).  It 
is  therefore  not  surprising  that  the  numerical  implementation  can  be  performed  in 
close  analogy  with  the  finite  element  method  using  efficient  numerical  tools  created 
within  the  last  decades. 

In  order  to  develop  this  analogy  formally,  the  basic  properties  of  the  finite  element 
method  will  be  briefly  outlined  at  this  point.  For  details  reference  is  made  to  the 
pertinent  literature  [2]. 

The  simple  static  element  in  Fig.  3a  is  chosen  as  a  representative  example.  It  has 
four  degrees  of  freedom,  here  the  four  independent  node  displacements  {u}m  = 
{’ti,U3,u3,U4}.  The  corresponding  node  forces  {p}m  =  {pi,pt,Ps,Pi}  are  then 
given  by 

{p}m  =  [*]m{«}m,  m  =  l,2...fiT,  (88) 

where  [A]m  is  the  local  element  stiffness  matrix  and  the  subscript  refers  to  the  actual 
element  number.  N  is  the  total  number  of  elements.  As  the  nodes  are  in  general 
common  to  more  elements,  it  is  here  convenient  to  introduce  a  global  degree-of- 
freedom  vector  {#}  defined  by 

{u}m  =  [L]m{U],  m  =  1,2...  fir,  (89) 

where  [L]m,  a  topology  or  connectivity  matrix  for  element  m,  is  an  extremely  sparse 
matrix  with  only  four  non-vanishing  elements  equal  to  unity.  The  mapping  de- 
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Fig.  3.  Analogy  to  finite  element  method:  (a)  finite  element  with  font 
degrees  of  freedom;  (b)  finite  wave  element  with  four  degrees  of  freedom. 


fined  by  Eq.  (89)  directly  reflects  the  connectivity  conditions  governing  the  globally 
assembled  element  system  in  question. 

At  this  point  a  variational  principle,  often  Hamilton’s  principle  of  stationary  energy, 
together  with  Eqs.  (S8),(89)  leads  to  the  following  linear  system  of  equations  having 
to  be  satisfied: 

WW  =  {*},  (90) 

where  [IT  j  is  the  global  stifihess  matrix 

(JT)  =  £>£l4.(L]m,  (81) 

m=l 

and  {R}  is  the  global  external  force  vector 

m  <«) 

mss  l 

The  local  vector  {r}m  represents  the  external  forces  acting  os  element  number  m, 
again  concentrated  in  the  nodes. 
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It  It  not  necessary  here  to  go  into  details  concerning  the  actual  implementation  of 
the  finite  element  equations,  but  it  it  important  to  note  the  structure  of  the  global 
system  of  equations,  Eq.  (90).  Since  the  topology  matrices  contain  only  seros  and 
ones,  the  matrix  multiplications  in  Eqs.  (90)  and  (91)  never  need  to  be  performed, 
but  can  be  replaced  by  a  set  of  indices  (pointers)  indicating  the  row  and  column  of 
the  global  stiffness  matrix  where  each  element  of  the  local  ones  have  to  be  added. 
The  actual  generation  of  the  global  stillness  matrix  can  therefore  be  performed  very 
efficiently  once  the  local  ones  have  been  determined. 

Returning  to  the  determination  of  the  depth- dependent  Green’s  function  in  stratified 
Media,  the  fact  that  the  horisontal  range  dependence  has  been  removed  by  the 
depth- separation  yields  the  possibility  of  effectively  representing  each  layer  as  a 
one  dimensional  ‘finite  wave  element’,  as  shown  in  Fig.  3b,  in  formal  analogy  to  the 
standard  finite  element  discussed  above.  The  degrees  of  freedom  for  this  element  are 
simply  the  amplitudes  >4~ ,  j4+ ,  1?~  and  27+  of  the  corical  waves  in  layer  number  m. 
These  are  conveniently  collected  in  a  local  degree-of- freedom  vector  (a(fc)}m,  where 
k  is  the  horisontal  wavenumber: 


{a(k)}m 


*.(*)  i 

AUk)  (  ’ 

*i(*>  J 


m  ^  1,2... IV. 


(93) 


If  the  kernels  for  the  field  parameters  involved  in  the  boundary  conditions  are  ex¬ 
pressed  in  vector  form  as 


{v(k,z)\ 


< 


w(i,z)  v 

**(M)  l 
*«(M)  f  ’ 
*r»(MMm 


m  =  1,2...JV, 


(94) 


the  following  matrix  relation  is  obtained  for  the  homogeneous  part  of  the  solution 
in  layer  number  m 

{»(*,  *)}m  =  [c(k,  i)]m{«(i)}m,  m  =  l,  2 . . .  N.  (95) 

The  local  coefficient  matrix  [c(i,  z)]m  is  a  function  of  the  horisontal  wavenumber  k 
and  the  depth  *.  Tot  the  homogeneous  layers  the  depth  is  involved  in  the  exponential 
functions  only,  and  in  these  cases  the  coefficient  matrix  can  be  factorised  as 

l<K*i*)]m  =  (d(*)]m[«(M))m,  (96) 

where  (d(k)]m  is  a  depth-independent  matrix  containing  only  simple  functions  of  k, 
and  [e(fr,  x]m  is  a  diagonal  matrix  containing  the  exponentials. 
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The  kernels  for  the  source  field  {v(Jfe,z)}m  ere  now  superimposed  with  the  kernels 
for  the  homogeneous  solution,  and  the  continuity  of  the  field  parameters  at  interface 
number  m  separating  the  layers  m  and  m  -f  1  can  then  be  expressed  as 

{v{k))Z  +  {v(k))Z  =  {v(k))Z+i  +  {v(k))Z+l,  m  =  1,2... N  -  1,  (97) 

where  the  depth  z  of  the  interface  has  been  removed  and  replaced  by  a  superscript 
indicating  the  interface  number.  If  Eq.  (97)  is  rewritten  as 

{«**)>£  -  M*)}£+1  =  {«(*)}£+ 1  -  W)}£,  m  =  1,2...  —  1,  (98) 

it  expresses  the  cancellation  of  the  discontinuity  in  the  source  field  by  the  discon¬ 
tinuity  in  the  homogeneous  solution.  The  interface  discontinuity  vector  (v(Jb)}m  is 
therefore  introduced  as 

w*)}m  =  {«(*)}£  -  m  =  1,2. .  .N  -  1,  (99) 

and  similarly  for  the  source  field  discontinuity  vector  {v(Jb)}m. 

In  order  to  assemble  the  local  equations,  Eq.  (98),  into  a  global  system,  the  globed 
degree- of- freedom  vector  {d(k)}  in  the  upgoing  and  downgoing  wavefield  amplitudes 
is  first  introduced,  defined  by  the  unique  local- to-global  mapping 

{a(*)}m  =  (S]m{A(*)},  m  =  1, 2 . . .  N.  (100) 


After  insertion  of  Eqs.  (95)  and  (100),  the  discontinuity  vector,  Eq.  (99),  takes  the 
form 

{«(*)}“  =  flcWClSl-  -  M*)S+ISU0  m  =  1, 2 . . . N  -  1.  (101) 


We  now  introduce  a  second  unique  mapping,  collecting  the  local  field  discontinuity 
vectors  {«(*)}"*  into  one  global  discontinuity  vector  (V(fc)}, 

{n*»=  xiprww.  (102) 

m=l 

which  after  insertion  of  Eq.  (101)  becomes 

<  vw>  =  £  trr  (w*)m»  -  w*)ctisi-+.)  {a<*>).  (um> 

mal 

Similarly  the  globed  source-field  discontinuity  vector  (V(ifc)}  is 

{?(*»  =  Em"  («*»s  -  <®w>:+.)  •  (104) 

m=l 
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Local  coefficient  matrix  [cj"  Pointer  matrix  [l]“ 


Global  coatflciant  matrix  [c] 


Fig.  4.  Mapping  between  local  and  global  coefficient 
matrices  by  mean*  of  row  and  colunn  pointers. 


The  global  cancellation  of  Eq.  (103)  by  Eq.  (104)  therefore  requires  the  following 
linear  system  of  equations  to  be  satisfied: 

IC(*)1M(A)}  =  -{?(*)},  (105) 

where  (C(fc)]  is  the  global  coefficient  matrix 

!<?(*)]  =  £pr  -  W*C+(S|m+.)  ■  (106) 

m=l 


As  can  be  observed,  the  global  system,  Eq.  (105),  is  very  similar  to  the  finite  element 
system,  Eq.  (90).  The  mapping  matrices  (T]m  and  [5]m  are  equivalent  to  the  topol¬ 
ogy  matrices  [£>]m  defined  for  the  finite  element  method  in  Eq.  (89).  However,  due 
to  the  fact  that  the  governing  boundary  conditions  for  the  wave  equation  are  not  set 
up  in  the  unknowns  directly,  but  in  derived  quantities,  two  sets  of  topology  matri¬ 
ces  are  needed  here:  one  for  the  local- to-global  mapping  of  the  degrees  of  freedom, 
i.e.  the  wavefield  amplitudes,  and  one  for  the  physical  parameters,  i.e.  displacement 
and  stresses  involved  in  the  boundary  conditions. 
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The  matrices  [S]m  and  [T]m  are  extremely  sparse,  containing  only  seros  and  ones. 
Since  the  mappings  of  Eqs.  (100)  and  (102)  are  unique,  the  corresponding  sum¬ 
mations  and  matrix  multiplications  in  Eqs.  (104)  and  (106)  need  never  actually  be 
performed  but  can  be  replaced  by  a  unique  set  of  pointers,  connecting  the  elements 
of  the  local  systems  with  those  of  the  global  system,  as  illustrated  in  Fig.  4.  As  is 
the  case  in  finite  element  programs,  the  topology  matrices  are  therefore  never  set  up 
in  the  actual  computer  code.  Their  formal  use  in  Eqs.  (104)  and  (106)  is,  however, 
very  convenient  in  the  general  fluid/solid/vacuum  case.  The  topology  matrix  [5]m 
is  set  up  to  include  only  the  non-vanishing  wavefield  amplitudes,  e.g.  in  the  upper 
and  lower  halfspaces  only  those  corresponding  to  upgoing  (A+,B+)  and  downgoing 
waves  respectively,  and  in  fluid  media  only  A ~  and  A+.  Similarly,  [T]m 
is  set  up  to  include  only  the  actual  number  of  boundary  conditions  at  interface  m, 
as  described  in  Sect.  3.5.  This  ensures  that  the  global  coefficient  matrix  is  a  square 
matrix  and  also  reduces  the  total  number  of  equations  and  unknowns  compared  to 
the  purely  solid  case. 

The  pointer  indices  defined  by  the  mappings  of  Eqs.  (100)  and  (102)  depend  solely 
on  the  number  of  unknowns  within  each  layer  and  the  boundary  conditions  at  each 
interface.  They  are  therefore  both  frequency  and  wavenumber  independent  and  can 
be  determined  a  priori. 

In  this  notation,  the  set-up  of  the  global  coefficient  matrix  requires  only  the  calcula¬ 
tion  of  the  elements  of  the  local  coefficient  matrices  [c(4)]m,  followed  by  the  indexed 
move  defined  by  the  mappings  given  in  Fig.  4.  The  subsequent  solution  of  Eq.  (105) 
then  yields  the  unknown  wavefield  amplitudes  in  all  layers  simultaneously. 

Although  the  global  system  of  equations,  Eq.  (105),  is  analytically  well  conditioned, 
apart  from  poles  corresponding  to  normal  modes  and  interface  waves,  its  numerical 
solution  is  not  necessarily  stable.  However,  by  choosing  the  mappings  such  that 
the  coefficient  matrix  becomes  block-diagonally  dominant  it  is  possible  to  ensure 
unconditionally  stable  solutions  using  gaussian  elimination  with  partial  pivoting. 
This  important  feature  will  be  discussed  in  the  following. 

The  difference  in  absolute  dimension  between  the  displacements  and  stresses  can 
y*eld  a  difference  of  several  orders  of  magnitude  between  the  coefficients  in  the  cor¬ 
responding  rows  for  both  the  local  and  global  systems.  As  is  well  known  in  such 
cases,  simple  gaussian  elimination,  with  or  without  partial  pivoting,  will  not  guaran¬ 
tee  unconditional  numerical  stability.  In  the  SAFARI  implementation,  the  coefficient 
matrices  are  therefore  scaled  to  make  all  elements  physically  dimensionless. 

Although  from  a  theoretical  point  of  view  the  origin  of  the  depth  axis  can  be  chosen 
arbitrarily,  its  choice  is  quite  critical  for  the  numerical  stability.  The  Airy  function 
solutions  for  the  inhomogeneous  fluid  layers  are  invariant  to  the  choice  of  origin, 
but  this  is  obviously  not  the  case  for  the  isovelocity  layers.  First  of  all,  a  common 
origin  for  all  layers  is  inconvenient  due  to  the  limited  allowed  arguments  to  the 
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exponential  functions.  A  local  origin  is  therefore  used  for  each  layer,  but  this  still 
does  not  ensure  numerical  stability  for  large  layer  thicknesses  and  large  horisontal 
wavenumbers  where  the  real  parts  of  arguments  to  tbe  exponential  functions  become 
significant  (the  evanescent  regime).  In  these  cases  argument  overflow  may  of  course 
appear,  but  even  before  this  happens  numerical  instability  may  appear  due  to  the 
fact  that  the  difference  between  the  growing  and  decaying  exponentials  exceeds  the 
machine  precision.  This  problem  can  be  avoided  by  factorising  out  the  exponential 
functions  growing  in  depth  [33],  but  unconditional  stability  can  also  be  obtained 
by  very  simple  means  not  requiring  any  additional  computational  effort,  just  by 
observing  the  physical  significance  of  the  problem  and  adopting  a  numerical  'trick' 
from  the  finite  element  technique.  This  is  the  key  to  the  efficiency  of  the  SAFARI 
code  and  will  therefore  be  discussed  at  this  point. 

The  physical  significance  of  the  evanescent  regime  for  a  certain  layer  is  that  the 
wavefield  in  the  absence  of  a  source  is  a  superposition  of  two  exponentially  decaying 
fields,  one  arising  from  a  non- vanishing  field  at  the  interface  above,  the  other  from 
the  interface  below.  Now  assume  that  sources  are  only  present  above  the  actual 
layer.  Then  the  field  decaying  away  from  the  lower  interface  (eai)  can  only  be  due 
to  'reflection'  from  this  interface  of  the  field  decaying  away  from  the  upper  interface 
(e~a*).  If,  however,  the  layer  is  very  thick  or  the  horisontal  wavenumber  is  large, 
the  amplitude  at  the  lower  interface  of  the  downwards  decaying  field  is  insignificant, 
and  the  upwards  decaying  field  will  be  practically  non-existent.  In  other  words,  the 
layer  containing  the  evanescent  field  will  behave  exactly  like  an  infinite  halfspace. 

It  is  obvious  that  numerical  stability  can  be  obtained  in  these  cases  if  th<  environ¬ 
mental  model  is  truncated  by  replacing  the  thick  layer  by  an  infinite  halfspace.  This 
would,  however,  introduce  significant  book-keeping  problems  and  rearrangements  of 
the  global  system  of  equations,  seriously  affecting  the  computation  time. 

Instead  it  is  possible  to  take  advantage  of  the  knowledge  about  the  vanishing  up¬ 
wards  decaying  wavefield  in  smother  way,  simply  adopting  a  numerical  trick  from 
implementations  of  the  finite  element  technique.  If  a  certain  node  displacement  is 
known  to  vanish  it  is  not  removed  from  the  degrees  of  freedom;  instead  the  corre¬ 
sponding  diagonal  element  in  the  global  stiffness  matrix  is  set  to  a  very  large  number. 
This  automatically  ensures  that  the  node  displacement  obtained  by  the  solution  will 
vanish  numerically. 

In  the  present  case  we  know  that  the  amplitude  of  the  upward-decaying  field  must 
vanish.  This  could  be  accomplished  by  placing  a  large  number  in  the  correspond¬ 
ing  diagonal  element  in  the  global  coefficient  matrix.  Again,  however,  complicated 
book-keeping  is  involved.  The  same  result  can  be  obtained  almost  free  of  charge 
in  the  following  way.  First  the  interface  above  the  layer  is  selected  as  the  local 
origin  for  the  depth  axis.  This  removes  the  exponential  functions  from  the  locad 
cocAdent  matrix  at  this  interface.  Then  the  argument  to  the  exponential  func¬ 
tion  is  truncated  at  a  value  close  to  the  maximum  allowed  for  the  actual  computer. 
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This  avoid*  overflow  problems  and  at  the  same  time  yields  very  large  numbers  in  the 
columns  of  the  global  coefficient  matrix  corresponding  to  the  upward-decaying  fields. 
It  can  now  be  shown  (20,21],  that  gaussian  elimination  with  pivoting  by  columns 
will  automatically  force  the  corresponding  amplitudes  to  vanish  numerically  if  the 
mappings  of  Eq.  (100)  and  (102)  are  chosen  such  that  the  global  coefficient  matrix 
becomes  block-diagonally  dominant.  In  essence  this  is  obtained  by  assembling  the 
local  systems  in  consecutive  order  and  placing  the  coefficients  involving  the  large  ex¬ 
ponential  functions  close  to  the  diagonal  of  the  global  matrix.  When  the  amplitude 
of  the  upward-decaying  field  is  thus  forced  to  vanish,  the  solution  obtained  for  all 
wavefield  amplitudes  in  the  other  layers  is  indistinguishable  from  the  one  obtained 
if  the  thick  layer  had  been  replaced  by  a  halfspace. 


Similar  arguments  can  be  used  to  show  that  also  in  the  case  of  sources  below  or 
within  the  thick  layer  the  SAFARI  technique  is  unconditionally  numerically  stable. 
Concerning  the  inhomogeneous  fluid  layers,  the  two  independent  solutions  Ai(()  and 
Bi(()  exponentially  decay  and  grow  respectively  for  £  — ►  oo.  The  stability  in  cases 
where  these  are  involved  is  therefore  ensured  in  exactly  the  same  way  as  described 
above. 


In  summary,  unconditionally  stable  solutions  are  obtained  simply  by  choosing  a 
proper  local  coordinate  system  and  a  proper  local-to-global  mapping  in  connection 
with  gaussian  elimination  with  partial  pivoting.  The  solution  technique  is  stable 
for  any  machine  precision.  Double  precision  is  therefore  not  required  for  stability 
purpose,  but  the  dynamic  range  is  of  course  limited  to  the  number  of  mantissa  digits 
available,  due  to  the  global  nature  of  the  solution.  Fbr  most  purposes,  however,  the 
6-7  orders  of  magnitude  offered  by  single  precision  is  more  than  enough. 


When  the  global  system  of  equations,  Eq.  (105),  has  been  solved,  the  kernels  in 
the  integral  representations  can  readily  be  evaluated  at  any  depth.  The  present 
technique  is  therefore  highly  efficient  in  cases  where  the  total  wavefield  is  to  be 
determined  at  many  depths,  e.g.  for  depth-range  contouring  of  single  frequency 
transmission  loss  and  for  synthetic  vertical  seismic  profiling  (VSP). 
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4.3.  RBPLBCTION  COBrPICIBNTS 

The  determination  of  plane- wave  reflection  coefficients  of  a  stratified  halfspace  is  of 
great  interest  in  both  underwater  acoustics  and  seismology  for  several  reasons.  It  is 
required  as  input  to  the  widespread  ray-tracing  codes  and  it  forms  the  basis  of  the 
synthetic  reflectivity  seismograms  commonly  used  in  crustal  seismology  [23].  Fur¬ 
ther,  the  plane-wave  reflection  coefficients  are  often  used  for  ocean  bottom- property 
inversion  schemes. 

The  full  wavefield  solution  technique  is  based  on  a  decomposition  of  the  total  wave- 
field  in  conical  waves  for  cylindrical  geometry  and  plane  waves  for  plane  cartesian 
geometry.  We  have  here  given  the  field  representations  in  cylindrical  geometry,  but 
the  similar  plane-geometry  representations  are  obtained  simply  by  replacing  the 
Hankel  transform  with  the  Fourier  transform;  the  kernels  are  identical  in  the  two 
cases  [21].  The  relation  between  the  horizontal  wavenumber  k  and  the  grazing  angle 
6  of  propagation  for  an  isovelocity  layer  is 

k  —  km  cos  6,  (107) 

where  the  medium  wavenumber  km  may  correspond  to  either  compressional  or  shear 
waves.  As  can  be  observed,  only  wavenumbers  k  <  km  correspond  to  real  angles. 


We  shall  now  demonstrate  that  the  plane-wave  r?lection  coefficient  Jt(9)  for  a  wave 
incident  from  above  at  angle  9  on  a  certain  interface  is  straightforwardly  found  by  the 
present  solution  technique.  The  layer  above  the  interface  is  replaced  by  an  infinite 
halfspace.  By  letting  the  source  field  be  a  plane  wave  incident  at  angle  9  and  with 
amplitude  A"(k)  =  A~(krn  cos  6)  at  the  interface,  the  solution  of  the  global  system 
of  equations  for  this  new  layered  problem  will  directly  yield  the  complex  amplitude 
of  the  reflected  plane  wave  A+(k)  =  A+{km  cos  5),  and  the  reflection  coefficient  can 
be  determined  by 


R(9)  = 


A+(km  cos  0) 
A~(km  cos  9)  ’ 


(108) 


Both  the  incident  and  reflected  waves  may  of  course  in  general  be  either  compres¬ 
sional  or  shear  waves.  As  the  global  solution  is  exact  to  within  machine  accuracy, 
the  same  is  the  case  for  the  reflection  coefficients,  i.e.  no  approximation  except  for 
the  environmental  model  is  involved. 


The  SAFARI-FIPR  module  computes  the  plane-wave  reflection  coefficients  for  inci¬ 
dent  compressional  waves  only,  since  a  standard  SAFARI  compressional  line  source  - 
placed  just  above  the  interface  -  is  used  to  produce  the  incident  plane  wave;  however 
the  calculated  reflection  coefficient  may  be  either  for  reflected  compressional  waves 
or  for  reflected  shear  waves. 


It  should  be  pointed  out  that  the  reflection  coefficients  calculated  in  this  way  are 
for  true  plane  waves.  Plane-wave  reflection  coefficients  are  difficult  to  determine 
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experimentally  [22],  which  should  be  kept  in  mind  when  comparing  synthetic  and 
experimental  reflection  data. 


4.3.  WAVENUMBER  INTEGRATION 


The  horisontal  wavenumbers  at  which  the  depth-dependent  Green's  functions  have 
to  be  determined  are  controlled  by  the  numerical  integration  scheme  used  for  evalu¬ 
ating  the  integral  representations  for  the  displacements  and  stresses.  These  Hankel 
transform  integrals  are  of  the  form 


G(r,z)=  J~  g{ktz)Jm(kr)kdkt 


(109) 


where  m  =  0  except  for  the  horisontal  displacement  and  shear  stress  where  m  =  1. 


lb  evaluate  these  integrals  numerically,  a  truncation  of  the  infinite  integration  in¬ 
terval  is  required.  Due  to  the  exponentially  decaying  kernels  for  k  oo  this  can  be 
done  to  any  degree  of  accuracy.  In  most  cases  a  truncation  point  10  to  20%  higher 
than  the  maximum  medium  wavenumber  in  the  problem  will  include  everything. 
There  are,  however,  exceptions  when  solid  layers  are  present.  In  these  cases  poles 
may  exist  due  to  evanescent  surface,  interface  or  plate  modes.  Especially  in  the  case 
of  thin  plates,  very  large  wavenumbers  are  involved.  It  is  therefore  important  to  be 
aware  of  the  behaviour  and  significance  of  these  waves  in  order  to  choose  a  proper 
integration  interval. 


Even  with  a  proper  truncation  there  sure  two  factors  complicating  the  numerical 
evaluation  of  the  Hankel  transforms.  One  is  the  kernel,  which  in  the  case  of  guided 
propagation  in  lossless  media  is  known  to  have  poles  on  the  real  axis  corresponding 
io  normal  modes  and  interface  waves  [6].  In  the  more  realistic  cases  where  vol¬ 
ume  attenuation  is  included,  these  poles  will  move  out  into  the  complex  plane,  in 
principle  making  real  axis  integration  possible,  but  the  kernel  will  in  most  cases  re¬ 
main  rapidly  varying.  The  second  complicating  factor  is  the  Bessel  function  Jm(kr), 
which  especially  for  long  ranges  r  oscillates  very  rapidly  in  k.  Therefore  an  accurate 
evaluation  of  the  Hankel  transform  will  usually  require  a  very  fine  sampling  in  the 
wavenumber  k. 


There  are  algorithms  available  for  directly  evaluating  the  Hankel  transform.  The  one 
developed  by  Tseng  [34],  however,  requires  the  evaluation  of  an  FFT  for  every  range 
selected,  but  even  more  importantly  it  requires  a  numerical  separation  parameter 
which  is  not  easily  selected.  The  so-called  Fast  Hankel  Transform  [35]  is  very  efficient 
for  relatively  smooth  kernels,  but  not  well  suited  for  the  rapidly  varying  kernels  of 
the  waveguide  problems. 

Nevertheless  it  has  been  shown  [15]  that  except  for  ranges  shorter  than  a  few  wave¬ 
lengths  and  very  steep  propagation  angles,  accurate  evaluation  can  be  obtained  by 
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an  asymptotic  approximation.  First  the  Bessel  function  is  expressed  in  terms  of  two 
Hankel  functions 

Jm(kr)  =  J  (*<i>t*r)  +  »<,”(*■•))  •  (HO) 

With  the  present  choice  of  time/frequency  transform  the  term  Hm\kr)  corresponds 
to  incoming  waves  which  are  only  important  for  representing  the  standing  wavefield 
at  very  short  ranges.  This  term  is  therefore  neglected,  and  Hm\kr)  is  replaced  by 
its  asymptotic  approximation 


fclim oH^(kr)=  yZLe-‘(*r-<m+*)H.  (in) 

The  inverse  Hankel  transform  then  takes  the  form 

G(r,z)=  J^‘i{m+knW  J“  9(k,*)'/k*-i'kdk.  (112) 

Due  to  the  desire  for  generality,  which  excludes  the  use  of  dedicated  quadrature 
schemes,  the  truncated  wavenumber  space  is  discretised  equidistantly  in  SAFARI  as 
follows: 

*i  =  *wi.  +  JA*.  f  =  0,l...(M- 1),  (113) 

where  M  is  the  total  number  of  sampling  points. 

The  vsdue  of  the  field  is  often  required  at  a  large  number  of  ranges  r,  in  particular 
in  connection  with  the  common  underwater  acoustics  problem  of  determining  trans¬ 
mission  loss  as  function  of  range.  In  these  cases  the  Fourier  integral  in  Eq.  (112) 
is  very  efficiently  evaluated  by  means  of  the  so-called  ‘Fast  Field*  approach  [15]. 
This  technique  is  therefore  applied  in  the  single-frequency  transmission-loss  pro¬ 
gram  SAFARI-FIP.  The  range  suds  r  is  discretised  as 


rj  =  rmi*  +  j Ar,  j  =  0, 1 . .  .(M  -  1),  (114) 

where  the  range  step  Ar  is  determined  by  the  relation 

ArAk  =  2w/M,  (115) 

and  M  is  an  integral  power  of  2.  The  following  discrete  approximation  of  Eq.  (112) 
is  then  obtained: 

G(ritx)»  V  [ff(*,,r)e-<r-«»,AfcV^]e-<(,^/M)t 

(116) 

where  the  summation  can  be  performed  by  means  of  an  FFT  yielding  the  field  at 
all  M  ranges,  Eq.  (114),  simultaneously. 
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It  is  wall  known  from  th«  use  of  the  FFT  for  time/frequency  transformations  that 
undersampling  in  one  domain  causes  aliasing  (wrap-around)  in  the  other  domain. 
This  is  also  the  case  here,  where  the  kernel  is  often  rapidly  varying  due  to  poles  and 
branch  points  close  to  the  integration  axis.  The  reason  is  (36)  that  the  evaluation  of 
Eq.  (116)  does  not  yield  G(r,  x)  but  instead  £nl-oo  (?(»•+»»£,  s),  where  R  =  2x/Ah, 
i.e.  a  sum  of  the  signals  in  all  range  windows  of  width  R  =  M  Ar .  As  can  be  observed, 
wrap-around  will  happen  from  both  sides  of  the  actual  interval  and  therefore  also 
from  ranges  smaller  than  rmlo.  The  negative  range  axis  does  not  give  any  significant 
wrap-around  due  to  the  fact  that  only  waves  travelling  in  the  positive  range  direction 
are  included.  However  if  rmiD  >  0  the  signals  in  the  interval  [0,  rmiB]  are  wrapped 
into  the  interval  [/t,rmin  +  A],  Consequently  the  maximum  useable  range  is  always 
R  =  M Ar,  independent  of  the  choice  of  rm|n. 

In  the  case  of  lossy  media,  the  aliasing  from  ranges  larger  than  rmi„  +  R  can  be 
reduced  by  choosing  R  so  large  that  the  signal  is  known  to  die  out  within  the  range 
window.  However,  we  see  from  Eq.  (115)  that  a  long  range  may  require  a  substantial 
number  of  wavenumber  samples,  especially  when  the  wavenumber  interval  is  fixed. 


The  aliasing  problem  can,  however,  be  solved  by  moving  the  integration  contour 
away  from  the  poles,  i.e.  by  smoothing  the  kernels.  According  to  Cauchy's  theorem 
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the  integrals  in  the  complex  plane  between  two  points  is  invariant  to  a  change  in  the 
integration  contour.  Therefore  Eq.  (112)  can  be  replaced  by 

G(r,z)  =  j  9{k^)^ke~irkdk,  (117) 

where  C  is  the  contour  shown  in  Fig.  5.  The  contour  consists  of  three  linear  sections 
Ci,  Ci  and  Cj,  where  the  vertical  sections  of  length  €  are  chosen  at  the  points 
where  the  wavenumber  axis  would  in  any  case  be  truncated.  If  these  points  are 
chosen  where  the  kernels  are  small,  i.e  ^(^in*xi  ^ ) V^mn  ^  0,  and 

f  <C  kmmx-km\n  then  the  contributions  from  the  vertical  sections  become  insignificant 
compared_to  the  integral  along  the  horizontal  section  defined  by  Jb  =  k  +  ie.  By 
inserting  k  in  Eq.  (117),  we  obtain 

G(r,  z)e~tr  ~  g(k  +  ie,  z)y/kTTee~irk  dk.  (118) 


This  integral  can  again  be  determined  by  means  of  an  FFT,  with  the  result 

oo 

^2  G(Tj,z)e-«r’+nR)  ~ 

n=-oo 

V'2*,ri 
M- 1 


x  [iK*!  +  *€>z)e  ira,iml£kky/ki  +  *ej  (H9) 


1=0 


or,  after  multiplication  with  eer> , 

G(rj,  z)  ~  »tr>  “i(fc-‘»rJ 

y/2irrj 


M-l 


X  |y(fcj  +  ic,  z)e  ‘p»i«,Afc  v/,tI  +  i(j  e  <(**07*0 
i=0 


-  Y'C(ri  +  n/Z,z)e 

n#0 


-tnR 
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It  is  now  clear  that  all  signals  wrapped  around  from  ranges  larger  than  rmjn  -f  R  will 
be  attenuated  by  at  least  e~tR.  On  the  other  hand  signals  wrapped  around  from 
ranges  smaller  than  rmln  will  be  amplified  by  at  least  etR.  As  was  the  case  for  the 
real  axis  integration  the  maximum  range  is  therefore  rm.x  =  R  also  for  the  offset 
contour  integration. 
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The  actual  value  of  e  is  not  extremely  critical.  However,  if  it  is  chosen  too  large 
the  contributions  from  the  two  vertical  parts  of  the  contour  may  become  significant. 
On  the  other  hand  a  too  small  value  will  require  a  very  large  number  of  sampling 
points.  For  most  practical  purposes  an  attenuation  of  the  wrap-around  by  60  dB  is 
more  than  enough.  The  corresponding  value  of  the  contour  offset  is 


3 

■Rloge 


2x(M  -  l)loge 


(^m*x  fcmln) 


(121) 


which  even  for  a  relatively  small  number  of  sampling  points  M  will  ensure  that 
e  <  kmmx  -  km in,  and  thus  yield  insignificant  contributions  from  the  vertical  con¬ 
tours. 


Whereas  the  Fast  Field  integration  technique  is  highly  efficient  for  single-frequency 
transmission-loss  calculations,  its  use  is  inconvenient  in  the  case  of  wideband  pulse 
calculations.  If  the  pulse  response  is  required  at  more  than  a  single  range,  the 
wavenumber  sampling  distance  AJb  would  have  to  be  frequency  independent  in  order 
to  satisfy  Eq.  (115).  Furthermore,  the  pulse  response  will  usually  only  be  required 
for  a  relatively  small  number  of  ranges,  so  that  direct  numerical  integration  for 
each  individual  range  is  feasible.  This  is  therefore  the  approach  taken  in  the  pulse 
program  SAFARI-FIPP. 

There  is  an  optional  choice  of  two  different  quadrature  schemes.  The  first,  which 
is  the  default  scheme,  performs  a  standard  trapezoidal  rule  evaluation  of  the  inte¬ 
gral  in  Eq.  (112),  still  with  equidistant  wavenumber  sampling.  The  trapezoidal  rule 
integration  approximates  the  integrand  by  a  function  varying  linearly  between  the 
sampling  points.  It  is  therefore  obvious  that  this  technique  is  only  applicable  out 
to  ranges  where  the  product  of  the  kernel  and  the  exponential  function  is  well  rep¬ 
resented  by  a  linear  function.  The  kernel  can  be  smoothed  by  moving  the  contour 
out  into  the  complex  plane  as  described  above,  but  the  exponential  fimction  varies 
rapidly  at  long  ranges.  To  ensure  that  the  exponential  function  alone  is  well  rep¬ 
resented  by  a  linear  function,  the  wavenumber  sampling  must  satisfy  the  following 
inequality  [36] 

A kR  <  (122) 

which  by  comparison  with  Eq.  (115)  translates  into  a  maximum  range  which  is 
much  shorter  than  the  one  obtained  by  the  Fast  Field  technique.  Here  it  should  be 
pointed  out,  however,  that  the  FFT  technique  has  a  degraded  accuracy  at  longer 
ranges,  which  is  insignificant  on  the  logarithmic  scale  used  for  representation  of 
transmission  loss  but  which  definitely  is  important  in  connection  with  wideband 
pulse  calculations.  For  accurate  pulse  calculations  the  maximum  ranges  for  the  two 
techniques  are  almost  identical  and  determined  by  Eq.  (122). 

It  is  possible,  however,  to  obtain  accurate  solutions  at  much  larger  ranges  by  applying 
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the  generalised  Filon  integration  scheme  [37]: 

f  _Afc_ 

SAp 
AJb 


/>>' 


W)dk  = 


[A  •  A^°' 
[  “  +  /(t)«s«‘l]  ,  A,  =  0, 


(123) 


where  A /  =  f(b)  -  f(a)  and  similarly  for  the  other  functions.  This  quadrature 
scheme  is  exact  for  linear  variations  of  /(Jb)  and  p(Jfe).  In  the  present  case  S  =  ~ir 
and  g(k)  =  k,  i.e.  a  linear  function  of  k.  The  Filon  scheme  is  therefore  applicable 
to  any  range  provided  that  the  kernels  are  well  represented  by  linear  interpolation 
between  the  sample  points.  On  the  other  hand  the  error  due  to  the  ‘nonlinearity’  of 
the  kernel  will  increase  with  range,  and  the  Filon  scheme  therefore  also  has  a  finite 
range  of  applicability.  Since  this  range  of  applicability  depends  on  the  smoothness 
of  the  kernel,  it  is  not  possible  to  give  any  specific  value.  Mallick  and  Fraser  [36] 
have  found  that  whereas  the  wavenumber  sampling  required  for  the  trapesoidal  rule 
of  integration  is  inversely  proportional  to  u >r,  the  Filon  scheme  requires  a  sampling 
which  is  approximately  inversely  proportional  to  y/uf.  The  additional  computations 
involved  in  the  Filon  scheme  are  insignificant  compared  to  the  total  calculation  time 
involving  the  determination  of  the  depth- dependence  of  the  field.  It  is  clear  that  this 
scheme  has  significant  computational  advantages  due  to  the  reduction  in  the  required 
sampling,  in  particular  if  combined  with  the  kernel-smoothing  contour  offset. 


For  all  the  numerical  integration  schemes  described  above  the  kernels  must  vary 
smoothly  between  the  sampling  points.  One  of  the  effects  of  a  too  rapidly  vary¬ 
ing  kernel  is  the  wrap-around,  but  as  demonstrated,  this  problem  can  be  overcome 
by  choosing  a  complex  integration  contour,  which  in  effect  smooths  the  kernels  by 
moving  the  contour  away  from  the  poles  and  branch  points.  If,  however,  the  in¬ 
tegration  interval  is  truncated  at  points  where  the  kernel  amplitude  is  significant, 
the  discontinuities  so  introduced  will  also  give  rise  to  wrap-around.  In  broadband 
calculations,  in  particular,  this  will  cause  artificial  arrivals  which  will  often  be  indis¬ 
tinguishable  from  true  arrivals.  These  discontinuities  therefore  must  be  smoothed. 
In  SAFARI-FIP  and  SAFARI-FIPP  this  is  done  by  extrapolating  the  kernels  by  a  first- 
order  Hermite  polynomial,  where  the  coefficients  are  chosen  such  that  the  kernel  and 
its  first  derivative  become  continuous  at  the  truncation  points. 


In  summary,  it  should  be  stressed  that  the  evaluation  of  the  wavenumber  integrals 
is  by  far  the  most  critical  part  of  the  full-wavefield  solution  technique.  Whereas  the 
depth-dependent  Green’s  function  is  exact  to  within  the  computer  accuracy  once  the 
environmental  model  has  been  chosen,  the  integration  requires  both  truncation  and 
discretisation  of  the  wavenumber  interval,  a  procedure  which  is  not  easily  automated 
since  the  significance  of  the  errors  introduced  is  highly  dependent  on  the  character¬ 
istics  of  the  actual  problem.  As  described  above  there  are  remedies  available  for 
reducing  the  errors  to  insignificance,  but  these  tools  have  to  be  selected  carefully, 
involving  a  significant  portion  of  physical  knowledge  and  experience.  The  best  rule 
of  thumb  that  can  be  given  for  the  use  of  this  type  of  model  is  therefore: 
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Never  be  satisfied  with  the  first  result  -  check  the  convergence! 


4.4.  FREQUENCY  INTEGRATION 

The  last  step  involved  in  determining  the  full  pulse  response  is  the  evaluation  of  the 
inverse  Fourier  transform,  Eq.  (8).  This  integral  has  to  be  evaluated  for  all  field 
parameters,  ranges,  and  depths  of  interest,  i.e. 


F(r,  zft)=  f  G(r,z,u>)etutdLj. 
J  —  OO 


(124) 


Due  to  the  fact  that  the  frequency  interval  is  predetermined  by  the  source  spectrum, 
i.e.  the  frequency  dependence  of  Sw,  the  frequency  axis  is  easily  truncated  and  the 
Fourier  integral  is  therefore  most  conveniently  evaluated  by  means  of  an  FFT,  i.e. 


F(r,z,t)~  Awe<a''"1»t*  £  [G(r,  z,u>i)eit"l‘,*u]  e«uWN\ 


(126) 


where  the  frequency  and  time  axes  have  been  discretized  as  follows, 

=  tmin  4"  j  At,  Ij  =  0, 1  .  .  .  (N  —  1), 

=  Wmin  +  /Aw,  /  =  0,  1  .  .  .  (N  -  1), 

and  the  samplings  satisfy  the  relation 


(126) 

(127) 


At  Aw  =  2  v/N. 


(128) 


Due  to  the  symmetry  relation  G(r,  z,-w)  =  G(r,z, w)  it  is  clear  that  the  time 
function  F(r,z,t)  determined  by  Eq.  (124)  is  real.  In  order  for  the  discrete  form, 
Eq.  (125),  to  yield  a  real  time  series  it  is  required  that  the  frequency  discretization, 
Eq.  (127),  be  symmetric  around  w  =  0.  There  sue  algorithms  available,  the  so-called 
real  FFT’s  (RFFT),  which  take  advantage  of  the  symmetry  relation  of  the  kernel 
and  therefore  only  require  input  of  G(r,  z,w)  at  the  N/ 2  discrete  frequencies 

u>t  =  /Aw,  f  =  0,l...tf/2-  1,  (129) 

and  still  yield  the  time  series  at  all  N  times,  F<q.  (126).  Due  to  the  computational 
savings,  the  RFFT  is  always  used  to  evaluate  the  time  series  in  SAFARI-FIPP. 
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The  discretisation  of  the  frequency  axis  again  introduces  WTap-around  due  to  the 
limited  time  window  T,  which  is  related  to  the  angular  frequency  sampling  Au>  by 

T  =  N&t  =  2t/Au;.  (130) 

As  was  the  case  for  the  wavenumber  integration,  the  wrap-around  can  be  eliminated 
-  or  at  least  reduced  -  in  two  ways.  One  way  is  to  smooth  the  kernel  G(r,z,u>)  by 
moving  the  integration  contour,  which  is  the  technique  called  complex  frequency  in¬ 
tegration  [36].  The  other  way  is  to  choose  T  so  large  that  the  whole  signal  is  included 
in  the  time  window.  Although  the  first  approach  definitely  has  computational  ad¬ 
vantages,  the  latter  is  the  one  chosen  in  SAFARI-FIPP.  Thir  is  due  to  the  transient 
nature  of  wideband  signals  which  makes  it  much  easier  to  perform  a  truncation  of 
the  time  window  than  of  the  range  window.  No  signal  can  arrive  to  a  receiver  at 
a  higher  speed  than  the  fastest  wave  speed  in  the  problem.  This  therefore  deter¬ 
mines  the  beginning  of  the  time  window  for  which  no  wrap-around  will  happen  from 
earlier  times.  The  length  T  of  the  time  window  is  more  difficult  to  determine,  but 
with  some  experience  in  conjunction  with  convergence  tests  it  becomes  relatively 
straightforward.  One  approach  is  to  make  an  initial  calculation  with  a  very  long 
time  window  but  with  a  narrow  source  spectrum.  This  will  of  course  yield  relatively 
wide  pulses,  but  the  overall  arrival  structure  can  be  determined;  the  time  window 
can  then  be  reduced  and  the  full  source  spectrum  included. 

In  SAFARI- rIPP  it  is  possible  to  calculate  the  pulse  responses  for  many  receiver 
ranges  and  depths  through  a  single  solution  pass,  and  a  fixed  stmt  of  the  time 
window  is  therefore  inconvenient.  Instead  one  can  choose  a  ‘running’  window  defined 

by 

tmin^)  =  r/c  P,  (131) 

where  r  is  the  range  and  cr  is  a  constant  called  the  reduction  velocity.  This  constant 
has  to  be  large  enough  that  <min(r)  is  smaller  than  the  arrival  time  of  the  fastest 
signal  in  order  to  avoid  wrap-around.  In  most  applications  cr  is  conveniently  chosen 
to  be  equal  to  the  fastest  wave  speed  in  the  problem. 
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5.  Installation  of  SAFARI 

The  SAFARI  package  consists  of  the  three  main  programs  SAFARI-FIPR,  SAFARI-FIP 
and  SAFARI-FIPP  and  a  number  of  subroutines,  most  of  which  are  common  for  the 
three  main  programs. 

All  output  from  the  SAFARI  modules  is  produced  in  graphics  form,  either  as  standard 
curve  pic;;  or  three-dimensional  contour  plots.  The  plotting  is  performed  by  two 
main  programs,  FIPPLOT  and  CONTUR,  which  have  to  be  called  after  execution  of 
the  calculating  programs. 

At  SACLANTCEN  the  programs  are  implemented  on  an  FPS164  Attached  Processor 
linked  to  a  VAX  8600.  The  programs  are  therefore  highly  vectorised,  and  a  number  of 
FPS164  library  subroutines  are  used.  Apart  from  these  and  a  few  critical  subroutines 
used  for  generating  the  kernel  of  the  Green's  function  written  in  APAL64  assembler 
plus  special  subroutines  for  asynchronous  I/O,  the  codes  are  written  entirely  in 
FORTRAN  77.  A  special  FORTRAN  77  library  with  simulations  of  these  special 
subroutines  is  available,  however.  The  codes  for  the  FPS164  and  the  VAX  are 
therefore  basically  identical.  The  fact  that  standard  FORTRAN  77  has  been  used 
makes  installation  possible  on  any  computer  satisfying  the  hardware  requirements. 


S.l.  PRECISION  REQUIREMENTS 

The  solution  algorithm  is  stable  for  any  machine  precision,  but  the  argument  of  the 
exponential  function  has  to  be  truncated  accordingly,  in  order  to  avoid  arithmetic 
overflow.  The  maximum  argument  is  set  in  each  of  the  main  programs  through  the 
variable  AM.  For  VAX  single  precision  a  value  of  65  is  used,  whereas  300  is  used  for 
64-bit  machines  such  as  the  FPS164  and  the  CRAY. 

Although  the  actual  precision  has  no  influence  on  the  stability,  it  affects  the  dynamic 
range  of  the  results,  but  double  precision  will  only  be  required  in  very  special  cases. 
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5.3.  MEMORY  REQUIREMENTS 


The  codes  will  run  in  their  default  version  on  any  machine  with  more  than  250  k 
words  of  memory,  where  one  word  corresponds  to  a  real  number  in  the  actual  pre¬ 
cision.  It  is,  however,  possible  to  reduce  the  array  sixes  to  allow  installation  on 
machines  with  less  memory.  The  array  sixes  are  basically  controlled  by  three  pa¬ 
rameters  defined  in  most  subroutines.  If  changes  are  made,  they  therefore  must  be 
made  in  all  source  files.  The  parameters  are: 

NLA :  Maximum  number  of  layers  in  environmental  model,  including  the  upper  and 
lower  halfspaces.  The  default  is  NLA  =  24. 

NP:  Maximum  number  of  wavenumber  samples.  In  SAFARI-FIPP  NP  is  also  the 
maximum  allowed  number  of  frequency  samples.  Due  to  the  use  of  an  RFFT 
for  the  frequency/time  transform,  the  corresponding  maximum  number  of 
time  samples  is  2NP.  The  default  is  NP  =  4096. 

NED:  Maximum  number  of  source  and  receiver  depths.  The  default  is  NRD  =  1001. 

It  is  obvious  that  the  array  sixes  may  also  be  increased  by  any  amount  allowed  by 
the  available  hardware. 


5.3.  SOURCE  FILES 


The  main  programs  are  placed  in  separate  source  files,  but  the  subroutines  have 
been  grouped  together  in  a  few  relatively  large  files  according  to  their  use.  In  the 
following  the  source  files  and  their  contents  are  described. 


SJEFIPR30.FOR  The  reflection  coefficient  program  SAFARI-FIPR  version  3.0, 
which  uses  the  SAFARI  solution  technique  to  determine  the  reflection  coefficients 
of  any  horixontally  stratified  fluid/solid  halfspace.  The  implemented  output  options 
are 

•  Modulus  and  phase  as  function  of  grazing  angle  for  a  selected  number  of 
frequencies. 

•  Modulus  and  phase  as  function  of  frequency  for  a  selected  number  of  grasing 
angles. 

•  Modulus  contoured  vs  frequency  and  grasing  angle. 

•  There  is  a  choice  between  calculation  of  the  P-P  reflection  coefficient  and  the 
P-SV  reflection  coefficient.  The  last,  of  course,  only  have  meaning  for  a  solid 
upper  halfspace. 
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SJEFIP30.FOR  The  trams  mission  loss  program  SAFARI- FIP  version  3.0,  which  for 
a  single  frequency  solves  the  wave  .equation  for  a  horisontally  stratified  environment. 
The  output  is  controlled  by  a  number  of  options.  Some  of  the  features  are 

•  Solid,  liquid  or  vacuum  isovelocity  layers  and  half  spaces.  Variable  velocity 
allowed  in  fluid  layers. 

•  Interface  roughness  (all  interface  types). 

•  Cylindrical  or  plane  geometry. 

•  Multiple  sources  (vertical  line  arrays). 

•  Multiple  receivers. 

•  Calculation  of  normal  stress  (pressure  in  fluids). 

•  Calculation  of  vertical  particle  velocities. 

•  Calculation  of  horisontal  particle  velocities. 

•  Depth  averaging  of  calculated  field  properties. 

•  Range/depth  contouring  of  field  properties. 

•  Complex  integration  contour  for  use  in  lossless  cases. 


SJEFIPP30.FOR  The  pulse  program  SAFARI-FIPP  version  3.0,  which  solves  the 
wave  equation  for  a  horizontally  stratified  environment  at  a  selected  number  of 
frequencies.  Some  of  the  features  ore 

•  Solid,  liquid  or  vacuum  isovelocity  layers  and  half  spaces.  Variable  velocity 
allowed  in  fluid  layers. 

^terface  roughness. 

•  Cylindrical  or  plane  geometry. 

•  Multiple  sources  (vertical  line  arrays). 

•  Multiple  receivers. 

•  Several  choices  of  source  pulse  shapes. 

•  Synthetic  hydrophone  signals. 

•  Synthe*5*  'ismograms.  These  can  be  plotted  separately  for  each  receiver  or 

.«*cke-  a  stacking  can  be  performed  in  either  range  or  depth. 

•  Dispersion  curves  for  normal  modes. 

•  Dispersion  curves  for  seismic  interface  waves. 

•  Complex  ration  contour. 
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FIPPLOT.FOR  The  main  program  FIPPLOT,  which  performs  ail  plotting  of  results 
except  for  contour  plots.  FIPPLOT  is  based  on  the  DISSPLA  plotting  system.  The 
plot  data  are  transferee!  through  two  files.  One  file  contain  titles,  labels,  axir  limits 
etc.  This  file  is  n  formatted  file  and  the  default  values  can  be  changed  by  the  user. 
The  other  file  is  a  formatted  file  which  contains  the  data  for  plotting. 

CONTUR.FOR  The  main  program  and  associated  subroutines,  which  perform  the 
contour  plotting  optionally  using  a  DISSPLA  compatible  plotting  package  or  the 
raster  package  UNIF.AS.  The  two-file  interface  to  the  calculating  programs  is  also 
used  here. 


FIPPSJE30.FOR  High-level  subroutines  used  by  all  programs. 


FIPTSJE30.FOR  High-level  subroutines  used  by  SAFARI-FIP. 


FIPR5JE30.FOR  Subroutines  used  only  by  SAFARI-FIPR. 


FIP5SJE30.FOR  Lower-level  subroutines  used  by  all  programs,  i.e.  basically  all 
subroutines  involved  in  determining  the  depth-dependent  Green’s  functions. 


FIPUSJE30.FOR  Subroutines  used  by  all  programs  for  generation  of  printout  and 
plot  files. 

APMATHSIM.FOR  Subroutines  simulating  the  APMATH64  library  of  the  FPS164 
Attached  Processor.  Also  contains  FORTRAN  equivalents  to  those  subroutines, 
which  in  the  FPSl64-version  are  written  in  APAL64  assembler. 


5.4.  COMPILING  THE  SAFARI  CODES 

The  compilation  is  straightforward,  since  no  special  options  are  needed.  Thus  for 
VAX/VMS, 


$  FOE/LIST 
4  FOE/LIST 
$  FOE/LIST 
$  FOE/LIST 
t  FOE/LIST 
$  FOE/LIST 
I  FOE/LIST 
t  FOE/LIST 


SJ1FIP30 

SJIFIPP30 

SJEFXPE30 

FIPPLOT 

COXTUE 

FIPPSJE30 

FIPSSJE30 

FIPESJE30 
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$  FOR/LIST  FIPTSJI30 
$  FOR/LIST  FIP9SJE30 
I  FOE/LIST  APNATHSIH 


Fbr  convenience  the  subroutines  are  placed  in  two  libraries: 


*  LI1R/CREATE  APNATH64  APNATHSIH 

I  LIB R/ CEE ATE  SAFLIB  FIPPSJE30,FIPSSJE30,FIPTSJE30,- 
FIP0SJE3O , FIFES JE30 


The  compilation  of  SAFARI  for  the  FPS164  is  described  in  Appendix  A. 


5.5.  LINKING  THE  SAFARI  CODES 

The  linking  of  the  five  executable  modules  is  performed  as  follows  under  VAX/VMS: 


$  LINK/EXE-FAXFIP30  SJKFIP30 , S AFLXB/LXfi , APH AT194/LIB 
f  LINI/EXK-TAXFIPP30  SJEFIPP30,SAFLIB/LIB,APNATH64/LIB 
t  LIVR/EXK-VAXFIPE30  S JEFIPE30 , SAFLIB/LIB , APNATHtf 4/LIB 
$  LINK  FIPPLOT , OISSPLA/LIB 
•  LINK  C0NT01, DISSPLA/LIB, ONIRAS/LIB 


where  DISSPLA.OLB  is  a  library  containing  the  DISSPLA  subroutines;  UNIRAS.OLB 
is  the  UNIRAS  library  if  available.  If  UNIRAS  is  not  available,  the  call  to  the  subrou¬ 
tine  MAINRAS  in  the  main  program  included  in  CONTUR.FOR  must  be  removed. 
DISSPLA  can  of  course  be  exchanged  with  any  other  plotting  system  available. 

The  linking  of  SAFARI  for  the  FPS164  is  described  in  Appendix  A. 


5.5.  RUNNING  THE  SAFARI  CODES 

The  programs  are  now  ready  to  run.  On  the  VAX  it  is  most  convenient  to  make  one 
command  file  with  parameter  input  for  each  program: 


FIPR.COM 


I  ASSIEH/BSER 
t  AS3I0N/9SER 
t  ASSIEN/OSER 
t  ASSI8N/USER 


'PI*. OAT  F0E001 
•Pl’.PLP  F0E019 
»P1».PLT  F0R020 
* PI '.COE  F0R028 


!  INPOT  DATA  FILE 

(PLOT  PAEANETEE8  (FIPPLOT) 

I PLOT  DATA  FILE  (FIPPLOT) 

I  CONTONE  PLOT  PAIAKETEE8 
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9  ASSIGM/USER  ’Pl’.BDR  F0A029 
<  MX  TAXFIPR30 
9  ASSIGM/USER  ’Pl'.PLP  F0E019 
9  ASSIGM/USER  *P1».PLT  F0R020 
9  RUN  FIPPLOT 

9  ASSIGM/USER  ’Pl’.CDR  FOKOS5 
9  ASSIGN/USER  ’Pl'.BDE  F01017 
9  AGE  CONTUR 


(CONTOUR  DATA  FILE 
!  EXECUTE  SAFARI "-FIP1 
I  ASSIGN  PLOT  FILES  TO 
! LOGICAL  MAKES  USED  BT  FIPPLOT 
{EXECUTE  FIPPLOT 
{ASSIGN  CONTOUR  FILES  TO 
(LOGICAL  MAKES  USED  BT  COMTUR 
{EXECUTE  CONTUR 


FIP.COM _ 

9  ASSIGN/USER  ’PI ’.DAT  F0R001 
9  ASSIGN/USER  ’Pl’.PLP  F0R019 
9  ASSIGN/USER  ’Pl’.PLT  F0R020 
9  ASSIGN/USER  ’Pl’.CDR  F0RO28 
9  ASSIGN/USER  ’Pl’.BDR  F0R029 
9  RUN  TAXFIP30 

9  ASSIGN/USER  ’Pl’.PLP  F0R019 
9  ASSIGN/USER  ’Pl’.PLT  F0R020 
9  RUN  FIPPLOT 

9  ASSIGN/USER  ’Pl’.CDR  F0R05S 
9  ASSIGN/USER  ’Pl’.BDR  F0R017 
9  RUN  CONTUR 


{INPUT  DATA  FILE 
{PLOT  PARAMETERS  (FIPPLOT) 
{PLOT  DATA  FILE  (FIPPLOT) 
{CONTOUR  PLOT  PARAMETERS 
! CONTOUR  DATA  FILE 
{EXECUTE  SAFARI -FIP 
{ASSIGN  PLOT  PILES  TO 
(LOGICAL  NAMES  USED  BT  FIPPLOT 
(EXECUTE  FIPPLOT 
{ASSIGN  CONTOUR  FILES  TO 
(LOGICAL  NAMES  USED  BT  CONTUR 
{EXECUTE  CONTUR 


FIPP.COM _ 

9  ASSIGN/USER  ’PI ’.DAT  F0R001 
9  ASSIGN/USER  ’PI’  .HP  F0R019 
9  ASSIGN/USER  ’Pl’.PLT  FOR020 
9  ASSIGN/USER  ’Pl’.CDR  F0R028 
9  ASSIGN/USER  ’Pl’.BDR  F0R029 
9  RUN  TAXFIPP30 
9  ASSIGN/USBR  ’Pl’.PLP  F0R019 
9  ASSIGN/USER  ’Pl’.PLT  F0R020 
9  RUM  FIPPLOT 


(INPUT  DATA  FILE 
(PLOT  PARAMETERS  (FIPPLOT) 
(PLOT  DATA  FILE  (FIPPLOT) 
(CONTOUR  PLOT  PARAMETERS 
(CONTOUR  DATA  FILE 
(EXECUTE  SAFARI-FIPP 
{ASSIGN  PLOT  FILES  TO 
(LOGICAL  NAMES  USED  BT  FIPPLOT 
(EXECUTE  FIPPLOT 


The  command  files  shown  here  of  course  have  to  be  extended  with  the  definition 
of  the  proper  directories.  However,  this  is  installation  dependent  and  therefore  not 
included  here.  The  preparation  of  command  files  for  running  SAFARI  on  the  FPS164 
is  described  in  Appendix  A.  Details  on  running  FIPPLOT  and  CONTUR  are  given 
in  Appendixes  B  and  G  respectively  . 
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5  -  Institution  of  SAFARI 


Assume  that  a  file  INPUT.DAT  has  been  prepared  for  SAFARI-FIP.  The  user  then 
only  needs  to  give  the  command 


•  trip  INPUT 


and  the  program  will  run.  All  files  related  to  the  actual  run  will  have  the  same 
names  except  for  the  extensions.  Similarly  in  batch: 


t  SUBMIT  PIP/pmMIMPUT 


The  preparation  of  input  files  is  described  separately  for  the  different  SAFARI  mod¬ 
ules  in  the  following  sections. 
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0.  Running  SAFARI-FIPR 

The  SAFARI-FIPR  module  calculates  the  plane- wave  reflection  coefficients  for  an  ar¬ 
bitrarily  stratified  fluid/solid  halfspace.  The  layer  above  the  interface  of  interest 
should  be  replaced  by  an  infinite  halfspace  with  the  properties  existing  just  above 
the  interface.  The  code  will  then  automatically  place  a  compressional  line  source 
in  the  upper  halfspace  and,  for  each  selected  graxing  angle  and  frequency,  use  the 
global  matrix  technique  to  calculate  the  complex  amplitude  of  the  resulting  upgoing 
plane  wave  in  the  upper  halfspace.  The  complex  reflection  coefficient  then  follows  di¬ 
rectly  by  division  with  the  complex  amplitude  at  the  interface  of  the  incoming  wave 
produced  by  the  source.  The  complex  reflection  coefficient  Rc  is  most  conveniently 
displayed  as  modulus  |Uc|  and  phase  tan~1(Im[lic]/Re[/ic])-  In  underwater  acous¬ 
tics  it  is  common,  however,  to  represent  the  modulus  by  the  corresponding  reflection 
lota  defined  by 

Ran  =  -20log|Rcj 

and  this  is  the  parameter  displayed  in  the  graphic  output  produced  by  SAFARI-FIPR. 

The  results  are  presented  in  graphic  form  as  curve  plots  showing  the  dependence  of 
the  reflection  loss  and  phase  shift  as  function  of  either  grazing  angle  or  frequency. 
Further,  there  is  an  optional  choice  of  reflection  loss  contours  vs  grasing  angle  and 
frequency. 

All  inputs  to  SAFARI-FIPR  are  read  from  the  file  currently  assigned  to  the  logical  file 
FOROOl.  Before  running  the  program  the  user  has  to  assign  the  file  containing  the 
input  data  to  this  logical  name.  It  is  most  convenient  to  include  this  assignment  in 
a  general  command  file  also  assigning  file  names  to  the  logical  names  of  the  output 
files,  as  described  in  Subsect.  5.6. 

In  the  following  the  preparation  of  the  input  files  will  be  discussed  in  detail  followed 
by  an  outline  of  the  necessary  numerical  considerations. 
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S.l.  INPUT  FELBS  FOR  SAFARI-FIPR 


The  Input  data  are  structured  in  9  blocks.  The  first  5  blocks,  shown  in  Table  2, 
specify  the  title,  options,  environmental  parameters,  together  with  the  desired  gras- 
ing  angle  and  frequency  sampling.  The  last  4  blocks,  outlined  in  Table  3,  contain 
axis  specifications  for  the  graphical  output.  Some  of  these  blocks  should  always  be 
included  and  others  only  if  certain  options  have  been  specified.  The  single  blocks 
and  parameters  are  described  in  detail  in  the  following. 


TtbJe  2 

Parameters  of  SAFARI-FIPR  input  files:  Calculation  and  environmental  parameters 


Block 

Parameter 

Units 

Limits 

I 

TITLE: 

title  of  run 

- 

<  80  char. 

II 

optl  opt 2  ... 

output  options 

- 

<  40  char. 

III 

ML:  number  of  layers,  ind.  halfspaces 

- 

ML  >  2 

0 

depth  of  interface 

m 

- 

CC 

compressional  speed 

m/s 

CC  >  0 

CS 

shear  speed 

m/s 

- 

AC 

compressional  attenuation 

dB/A 

AC  >  0 

IS 

shear  attenuation 

dB/A 

AS  >  0 

10 

density 

g/cm* 

10  >  0 

16 

rms  value  of  interface  roughness 

m 

- 

CL:  correlation  length  of  roughness 

m 

CL  >  0 

IV 

rum 

minimum  frequency 

Hs 

FMIM  >  0 

FRIX 

maximum  frequency 

Hs 

FMAX  >  0 

N1F1 

number  of  frequencies 

- 

'•onum  >  i 

NFtJU 

plot  output  increment 

- 

N70V  >  0 

V 

AMU 

minimum  grasing  angle 

deg 

AMU  >  0 

AHAX 

maximum  grasing  angle 

deg 

AMAX  >  0 

MIAN 

number  of  grasing  angles 

- 

MIAN  >  1 

NAOO 

plot  output  increment 

NAOO  >  0 
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Table  3 

Parameters  of  SAFARI-FIPR  input  files:  Plot  parameters 


Block 

Parameter 

Units 

Limits 

VI1 

ALEF: 

left  border  of  angle  axis 

deg 

— 

ARIG: 

right  border  of  angle  axis 

deg 

- 

ALIN: 

length  of  angle  axis 

cm 

ALEN  >  0 

AINC: 

distance  between  tick  marks 

deg 

AINC  >  0 

HALO: 

lower  border  of  R-loss  axis 

dB 

— 

RAUP: 

upper  border  of  R-loss  axis 

dB 

- 

RFLN : 

length  of  loss  and  phase  axes 

cm 

RFLN  >  0 

RFIN: 

R-loss  axis  tick  mark  interval 

dB 

RFIN  >  0 

VII* 

FLEF: 

left  border  of  frequency  aria 

Hz 

FRIG: 

right  border  of  frequency  axis 

Hz 

- 

FLEN: 

length  of  frequency  axis 

cm 

FLEN  >  0 

FINC: 

distance  between  tick  marks 

Hz 

FINC  >  0 

RFLO: 

lower  border  of  R-loss  axis 

dB 

— 

RFUP: 

upper  border  of  R-loss  aria 

dB 

- 

RFLN: 

length  of  loss  and  phase  axes 

cm 

RFLN  >  0 

RFIN: 

R-loss  axis  tick  mark  interval 

dB 

RFIN  >  0 

VIII3 

ALEF: 

left  border  of  angle  axis 

deg 

ARIG: 

right  border  of  angle  axis 

deg 

- 

ALEN: 

length  of  angle  axis 

cm 

ALEN  >  0 

AINC: 

distance  between  tick  marks 

deg 

AINC  >  0 

FRLO: 

lower  border  of  frequency  axis 

Hz 

FRLO  >  0 

PROP: 

upper  border  of  frequency  axis 

Hz 

FRUP  >  0 

OCLN: 

length  of  one  octave 

cm 

OCLN  >  0 

NTKM: 

number  of  tick  marks  per  octave 

- 

NTKM  >  0 

ZHIN: 

minimum  contour  level 

dB 

— 

ZNAX: 

maximum  contour  level 

dB 

- 

ZINC: 

contour  level  increment 

dB 

ZINC  >  0 

IX4 

TLEF: 

wave  speed  at  left  border 

m/s 

_ 

TRIG: 

wave  speed  at  right  border 

m/s 

- 

TLEN: 

length  of  wave  speed  axis 

cm 

TLEN  >  0 

TINC: 

wave  speed  tick  mark  distance 

m/s 

TINC  >  0 

0T0P: 

depth  at  upper  border 

m 

- 

DTLO: 

deptit  at  lower  border 

m 

- 

DTLN: 

length  of  depth  axis 

cm 

DTLN  >  0 

DTIN : 

depth  axis  tick  mark  interval 

m 

DTIN  >  0 

1  Only  for  WOO  >  0.  1  Only  for  NAO0  >  0.  3  Only  for  option  C.  4  Only  for  option  Z. 
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mmmmm  -  options 


Mmt 

Required: 

$r*rt**: 

Fhrm*t: 

Dmtriptioa: 


A-/.." 

pptt  opt2  opt3  ... 

free.  Maximum  40  alphanumeric  characters. 

Scrips  of  alphanumeric  characters  specifying  the  desired 
options.  The  following  are  implemented  at  present: 


C :  Loss  contours  plotted  vs  frequency  and  gracing 
angle.  : 

lit  This  is  the  default  option,  which  therefore  never 
needs  to  be  specified.  It  calculates  the  P-P 
reflection  loss  as  Junction  of  angle  and  frequency. 

Ft  Phase  angle  of  reflection  coefficient  {dotted  in 
addition  to  the  de&nk  reflection  loss. 

8:  This  option  replaces  the  default  P-P  wave 
redaction  coefficient  by  the  F-SV  wave  reflection 
coefficient. 

Z:  Plot  of  velocity  profiles. 
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block  i  - 

TITLE 

Lines:  : 

l 

R/tguirtd: 

Always 

Sjfntax: 

TITLE 

Fbrmnt: 

A80 

Dmcription: 

TITLE:  Title  of  run.  Maximum  80  characters. 
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BLOCK  m  ~  ENVIRONMENTAL  DATA 

Lines:  >3 

Required:  Always 

‘  Syntax: 

ML 

D(l)  CC(1)  CS(1)  AC(l)  AS(1)  RD ( 1 )  RG(1)CL(1) 
D(2)  CC(2)  CSC2)  AC(2)  AS(2)  R0(2)  RG<2)  CL(2) 

D(NL)CC(NL)CS(NL)AC(HL)ASCNL)RO(NL)RG(NL)CL(HL) 


Description: 

NL:  Number  of  layers,  including  the  upper  and  lower 

halfspaces.  These  should  always  be  included,  even 
in  cases  where  they  are  vacuum.  In  this  regard  it 
should  be  noted  that  the  reflection  coefficient  is 
always  calculated  for  the  uppermost  interface, 
with  the  incoming  wave  being  a  downward- 
going  plane  P*wav«  In  the  Upper  halfspace.  This 
haHspace  is  therefore  not  alh^sred  tofe*  vacuum. 

DO :  Depth  r  an  m  of  upper  boundary  of  layer  or 

halfspace.  The  reference  depth  can  be  choosen 
arbitrarily,  and  DO  is  allowed  to  be  negative^ 

For  layer  no.  1,  i.e.  the  upper  halfspace,  this 
parameter  is  diunmy. 

CC():  Velocity  ce  ofcompressional  waves  in  m/s.  If 

specified  to  0.0,  the  layer  or  halfspace  is  vacuum. 

CSC):  Velocity  of  shear  waves  in  m/s.  In  order  to 
.  be  physically  meaningful,  it  is  required  that 
c,  <  VO.thce.  If  specified  to  0.0>  the  layer  or 
halfspace  is  fluid.  If  c,  <  0,  then  |c,J  represents 
the  compre8slonal  velocity  at  the  bottom  of  the 
:  actual  layer,  which  is  treated  as  fluid  with  l/c(r)* 

varying  linearly  with  depth. 

AC() :  Attenuation  of  compressional  waves  in  dB/A. 

If  the  layer  is  fluid,  and  AC()  is  specified  as  0.0, 
then  an  empirical  water  attenuation  is  used. 
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BLOCK  EC  -  (Continued) 

Description: 

A8():  Attenuation  jt  of  shear  wave*  in  dB/ A.  Fbr  the 
attenuations  to  be  physically  meaningful,  it  is 
required  that 


Density  p  in  g/cm3. 

RMS  roughness  of  into  face  in  m.  RG  (  J.  )  is 
dummy.  If  RG  >  0  the  KirchhoiF  approximation 
it  used  for  the  scattering  loss,  i.e.  CL  -+  oo. 

This  is  computationally  faster  than  doing  the 
full  scattering  theory,  but  it  is  only  applicable  for 
small  roughnesses  and  large  correlation  lengths. 
To  invoke  the  full  non-Kirchhoff  scattering  JIG 
must  be  specified  equal  to  the  negative  rms 
roughness,  i.e.  RG  <  6. 

Correlation  length  of  roughness  in  m  This  :Wm: 
parameter  is  only  required  if  nun*KirchhoiF 
scattering  has  been  selected,  i.e.  RG  <  0  ”  • 
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BLOCK  IV  - 

FREQUENCY  DATA 

Lines: 

1 

Required: 

Always 

?  Syntax: 

FHIS  FMAX  NRFR  HFOU 

Format: 

Jfree 

Description: 

FMIN :  Minimum  frequency  in  Hz. 

FHAX:  Maximum  frequency  in  Hz. 

NRFR:  Number  of  frequencies  (>  1).  If  the  contour 

option  C  was  specified  the  samples  will  be  placed 
equidistantly  on  a  logarithmic  frequency  axis; 
otherwise  the  sampling  will  be  equidistant  on  a 
linear  frequency  axis. 

HFOU:  IfNFOU  >  0  a  plot  of  reflection  loss  (surd 

phase  shift  if  option  P  was  chosen)  as  function 
of  grazing  angle  is  generated  for  each  HFOU 

frequency.  If  NFOU  =x  0  no  plots  as  function  of 
grazing  angle  will  be  produced, 
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BLOCK  V  ~ 

GRAmO  ANGLE  DATA 

lines; 

1 

Required: 

Always 

Syntax: 

AMIN  AMAX  NRAN  NADU 

Pbrmet:  - 

FVee 

Description: 

AMIN:  Minimum  grazing  angle  in  degrees. 

AMAX:  Maximum  grazing  angle  in  degrees. 

NRAN:  Number  of  angles. 

NAOU:  If  NAOU  >  0  a  plot  of  reflection  loss  (and  phase 
shift  if  option  P  was  chosen)  as  function  of 

frequency  is  produced  for  each  NAOU  angle.  If 

NAOU  =  0  no  plots  of  reflection  lose  as  function 
of  frequency  are  generated. 
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BLOCK  VI  - 

PLOT  DATA  /  REFLECTION  LOSS  AND  PHASE  VS 
GRAZING  ANGLE 

tinea:  ;  v 

2 

Required: 

If  NFOU  >  0 

::  Syntax: 

ALEF  ARIG  ALEN  AINC 

RALO  RAUP  RALN  RAIN 

Format:  T  • 
Description: 

Free 

ALEF :  Gracing  angle  in  degrees  corresponding  to  left 

border  of  plots. 

ARIG:  Right  border  of  grazing  angle  axis. 

ALEN:  Length  of  angle  axis  in  cm. 

AINC:  Distance  in  degrees  between  tick  vks  on 
grazing  angle  axis. 

RALO :  Reflection  loss  in  dB  corresponding  to  lower 

border  of  plot.  Fbr  phase  shift  plots  the  phase 
angle  axis  will  automatically  cover  the  interval 
-180°  to  180°.  This  interval  Can,  howevert  be 
changed  as  described  in  Appendix  B. 

RAUP;  Upper  border  in  dB. 

RALN :  Length  of  reflection  loss  and  phase  shift  axes  in 
cm. 

RAIN :  Distance  in  dB  between  tick  marks  for  loss  axis. 
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BLOCH  VB  ~ 

PLOT  DATA  /  REFLECTION  LOSS  AND  PHASE  VS 
FREQUENCY 

..  Lines:  ■ 

2 

Required: 

If  NAOU>  0 

Syntax: 

PLEF  PRIG  FLEN  FINC 

Fbrmnt; 

RFLO  RFUP  RFLN  RFIN 

Free 

Description: 

FLEF :  Frequency  in  H*  corresponding  to  left  border  of 
plots. 

FRIG:  Right  border  frequency  axis. 

/.FLEN :  Length  of  frequency  axis  in  cm. 

FINC:  Distance  in  He  between  tick  marks  on  frequency 
alia. 

RFLO :  Reflection  loss  in  dB  corresponding  to  lower 

border  of  plot .  For  phase  shift  plots  the  phase  ' 

angle  axis  will  automatically  cover  the  interval 
-180°  to  180°.  This  interval  can,  however,  be 
changed  as  described  in  Appendix  B. 

RFUT  i  Upper  border  of  reflection  loss  axis  in  dB. 

RFLN :  Length  of  reflection  loss  and  phase  shift  axes  in 
cm. 

RFIN :  Distance  in  dB  between  tick  marks  on  loss  axis. 
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blo  ck  van  - 

PLOT  DATA  /  ANGLE-FREQUENCY  CONTOURS 

Lines; 

3 

Requind: 

If  option  C  was  specified. 

Syntax: 

ALEF  AAIG  ALEN  AINC 

FRLO  FRDP  OCLN  NTKM 

ZNIN  ZMAX  ZINC 

Fbrm*t: 

Description: 

Free 

ALEF :  Gracing  angle  in  degrees  corresponding  to  left 
border  of  contour  plots. 

ARIG:  Right  border  of  grazing  angle  axis  on  contour 
plot. 

ALEN :  Length  of  grazing  angle  axis  in  cm. 

AINC :  Distance  in  degrees  between  tick  marks  on  angle 
axis. 

FRLO  :  Frequency  in  Hz  corresponding  to  lower  border  of 
plot.  Since  the  frequency  axis  will  be  logarithmic, 
FRLO  >0. 

FRUPi  Frequency  corresponding  to  upper  border  of  plot; 
OCLN :  Size  parameter  specifying  number  of  cm  per 
octave  on  frequency  axis, 

NTKH :  Number  of  tick  marks  per  octave. 

ZNIN:  Lower  limit  in  dB  of  reflection  loss  interval  for 
which  contour  lines  are  to  be  plotted. 

ZMAX:  Upper  limit  in  dB  of  reflection  loss  interval  for 
which  contour  lines  art  to  be  plotted. 

ZINC:  Distance  in  dB  between  contour  lines. 
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PLOT  DATA  f  VELOCITY  PROFILE 

-  Hum 

2 

If  option  Z  w&»  specified  in  Block  II. 

Sjmfx: 

VUP  ?MG  VWW  ¥IltC 

'  ^;i::  v::: -•  •  -  ■  :  • 

DVUP  0VLO  DVI.M  DVIM 

En* 

Dmwptioa: 

VLEF: 

Left  border  value  of  velocity  in  m/s. 

VRIG: 

Right  border  value  of  velocity  in  m/s. 

VLEN: 

Length  of  velocity  axis  in  cm. 

VI8C: 

Distance  between  tick  marks  in  m/s. 

D¥UP< 

Depth  at  upper  border  of  plot  in  an. 

v\>  "\>  .  ... 

DVLO  .* 

Depth  at  lower  border  of  plot  in  *». 

-  ..  .  .  . 

DVLN: 

Length  of  depth  axis  in  cm. 

x  ;'••••  ••  ■•••  "•••’ 

| ‘  . ' 

DTIM; 

Distance  in  m  between  tick  mirks  on  depth  axis. 
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9.2.  NUMERICAL  CONSIDERATIONS 

Since  the  determination  of  plane- wave  reflection  coefficients  does  not  involve  any  nu¬ 
merical  evaluation  of  integral  transforms,  the  only  numerical  problem  is  related  to 
the  discretisation  of  the  environmental  model.  When  modeling  real  data  where  only 
a  few  environmental  data  are  known,  it  is  obvious  that  the  model  has  to  be  re-run 
several  times  in  order  to  fit  a  discrete  stratified  model.  On  the  other  hand  the  accu¬ 
racy  requirements  are  very  limited,  and  there  are  no  purely  numeric  considerations 
involved.  This  is,  however,  the  case  when  the  SAFARI-FIPR  reflection  coefficients 
have  to  be  compared  to  those  obtained  by  other  numerical  models,  which  is  very 
common  in  relation  to  the  verification  of  new  algorithms.  Some  reflection  coefficient 
codes,  e.g.  those  based  on  finite  differences,  assume  a  linearly  varying  sound  veloc¬ 
ity  profile  between  the  sampling  points,  whereas  SAFARI  assumes  either  isovelocity 
layers  or  a  sound  speed  varying  non-linearly  between  the  sampling  points,  Eq.  (5). 
It  is  therefore  necessary  that  the  number  of  profile  sampling  points  is  so  large  that 
the  differences  between  the  interpolation  techniques  has  insignificant  effect  on  the 
calculated  sound  field.  It  is  not  possible  to  giv  '  any  specific  rules  in  this  regard  since 
the  effect  depends  on  the  frequencies  and  grazing  angles  of  interest.  If  isovelocity 
layers  are  used  to  represent  the  profile,  a  layer  thickness  of  J-A  will  usually  be  suf¬ 
ficient.  For  the  inhomogeneous  fluid  layers,  however,  the  thickness  can  be  chosen 
much  larger.  In  any  case  the  profile  discretization  should  be  controlled  by  checking 
the  convergence. 


«.s.  SAFARI-FIPR  EXAMPLES 

In  this  section  a  few  characteristic  examples  will  be  given  of  the  use  of  SAFARI-FIPR 
for  determining  plane-wave  reflection  coefficients.  The  examples  here  have  been 
chosen  from  underwater  acoustics,  but  problems  from  seismology  and  ultrasonics 
are  treated  in  an  identical  manner.  Several  other  examples  cm  be  found  in  [20-27]. 

The  environment  considered  in  the  present,  and  in  all  following  examples  for  the 
transmission  loss  and  pulse  models,  is  characterized  by  a  water  depth  of  100  m  and 
a  bottom  consisting  of  a  silt  layer  of  20  m  thickness  overlying  a  homogeneous  sub¬ 
bottom  of  coarse  sand.  The  material  properties  of  the  bottom  are  given  in  Table  4. 
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Table  4 

Material  parameters  for  ocean  bottom 


Layer 

Thickness 

Wave 

speed 

Attenuation 

Density 

P 

H 

Cc 

c. 

7c 

7. 

(m) 

(m/s) 

(m/s) 

(dB/A) 

(dB/A) 

(g/cn>S) 

sat 

20.0 

400 

0.2 

0.5 

1.8 

Sand 

oo 

600 

0.1 

0.2 

2.0 

■  6.3.1.  SAFARI-FIPR  case  1:  Reflection  coefficients  vs  grazing  angle 

This  test  case  concerns  the  determination  of  the  plane-wave  reflection  coefficient  of 
the  ocean  bottom  as  function  of  grazing  angle  at  a  single  frequency  of  50  Hz.  The 
data  file  is  given  in  Table  5. 

When  running  SAFARI-FIPR  with  the  above  data  file,  3  plots  are  produced  (Fig.  6). 
Figure  6a  shows  the  calculated  reflection  loss  with  the  velocity  profile  plot  inserted. 
Note  that  the  loss  is  very  small  for  grazing  angles  less  than  the  critical  angle  for  the 
sub-bottom,  9cJ  =  cos"1(1500/1800)  =  33.6°.  At  a  frequency  of  50  Hs  no  significant 
transition  happens  at  the  sediment  critical  angle,  0cl  =  cos-1  (1500/1600)  =  20.4°. 
The  main  effect  of  the  low  shear  speeds  is  the  non- vanishing  reflection  loss  at  grazing 
angles  smaller  than  critical.  Figure  6b  shows  the  associated  phase  shift  as  function 
of  grazing  angle.  Again  the  transition  at  the  sub-bottom  critical  angle  is  appar¬ 
ent.  Here  it  should  be  noted  that  no  phase  unwrapping  is  performed;  therefore  the 
discontinuous  behaviour  at  33.6°  in  fact  represents  a  360°  phase  change. 


■  6.3.2.  SAFARI-FIPR  case  2:  Loss  contours  vs  frequency  and  single 

In  this  example  we  calculate  the  reflection  loss  as  a  function  of  grazing  angle  over 
4  frequency  octaves  from  20  to  320  Hz  and  generate  a  contour  plot.  The  environment 
is  the  same  as  case  1,  and  the  data  file  is  given  in  Table  6. 

The  contour  plot  generated  by  SAFARI-FIPR  for  this  case  is  shown  in  Fig.  7a.  Fig¬ 
ure  7b  shows  the  plot  obtained  when  the  contour  program  CONTUR  is  re-run  with 
the  UNIRAS  option,  as  described  in  Appendix  C.  The  dark  areas  indicate  high  loss, 
the  light  ones  low  loss.  Note  that  the  contour  levels  have  been  changed  to  a  spacing 
of  2  dB  in  Fig.  7b. 

Note  that  at  the  low  frequencies  the  reflection  properties  are  dominated  by  the 
critical  angle  9e *  =  33.6°  of  the  sub-bottom,  whereas  for  the  higher  frequencies 
a  series  of  'resonance  ridges’  of  high  loss  appear,  asymptotically  approaching  the 
sediment  critical  angle,  9C\  =  20.4°.  This  behaviour  is  typical  for  layered  media. 
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Table  5 

SAFARI-FIPR  case  1  data  file 


Data  file 

Description 

Notes 

SiFAEl-FIPR  case  1 

P  Z 

title 

options:  P  -  phase  plot,  Z  -  velocity  profile 

1,2 

3 

0  1800  00010 

number  of  layers 

layer  1  (upper  halfspace,  water) 

3 

0  1800  400  0.2  0.8  1.8  0 

layer  2  (sediment) 

4 

20  1800  800  0.1  0.2  2.0  0 

layer  3  (sub-bottom) 

4 

80  80  1  1 

frequency  50  Hi,  NF00  =  1 

5 

0  90  181  0 

angle  interval  0  to  90°,  181  samples 

5,6 

0  90  20  30 

angle  axis  for  plots 

6 

0  16  12  6 

loss  axis  for  plots 

6 

0  2000  10  1000 
-20  40  10  20 

velocity  axis  for  profile  plot 
depth  axis  for  profile  plot 

1  In  order  to  obtain  plots  of  the  phase  shift  in  addition  to  the  default  plot  of  the  reflection 
loss,  option  P  has  been  specified. 

*  Option  Z  is  specified  to  get  a  plot  of  the  velocity  profile.  This  option  requires  that  the 
axis  specifications  are  given  at  the  end  of  the  data  file.  In  the  present  case  the  profile 
is  plotted  from  a  depth  20  m  above  the  bottom  to  40  m  below. 

3  Note  that  the  100  m  thick  water  column  has  been  replaced  by  an  infinite  halfspace 
with  the  properties  of  the  water.  For  convenience  the  ocean  bottom  has  been  selected 
as  origin  for  the  depth  axis.  This  is,  however,  arbitrary,  and  100  m  could  be  added  to 
all  depths  without  influencing  the  results. 

4  Since  both  interfaces  are  plane  the  rms  roughness  is  specified  as  0  in  all  cases,  and  the 
correlation  lengths  have  been  omitted. 

4  By  specifying  HF00  =  1,  plots  of  the  reflection  data  as  a  function  of  grasing  angle  will 
be  produced  for  each  frequency,  in  this  case  of  course  only  one.  NiOO  — .  0  indicates 
that  no  plots  of  reflection  loss  as  a  function  of  frequency  are  requested. 

*  The  reflection  coefficient  will  be  calculated  in  181  equidistantiy-spaced  points  in  the 
angular  interval  0  to  90*,  and  the  plots  will  show  the  reflection  loss  and  phase  shift 
over  the  whole  interval.  The  loss  axis  is  specified  to  cover  the  interval  0  to  15  dB. 
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(a)  REFLECTION  COEFFICIENT 


(b)  REFLECTION  COEFFICIENT 


Fig.  6.  Results  for  SAFARI-FIPR  case  1:  (a)  reflec¬ 
tion  loss  as  function  of  grasing  angle,  with  velocity 
profile  inserted;  (b)  phase  shift  as  function  of  grac¬ 
ing  angle 
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Table  6 

SAFARI-FIPR  case  2  data  file 


Data  file 

Description 

Notes 

SAFARI-FIPR  case  2 

C 

title 

option  C  -  contour  plot  of  reflection  loss 

1 

3 

0  1500  00010 

0  1600  400  0.2  0.5  1.8  0 

20  1800  600  0.1  0.2  2.0  0 

20  320  41  0 

number  of  layers 

layer  1  (upper  halfspace,  water) 

layer  2  (sediment) 

layer  3  (sub-bottom) 

frequency  interval  20  to  320  Hz,  41  samples 

2 

0  90  181  0 

angle  interval  0  to  90°,  181  samples 

2 

0  90  20  30 

angle  axis  for  contour  plot 

3 

20  320  3  1 

frequency  axis  fo’  contour  plot 

3 

1  10  1 

contour  levels 

4 

1  Option  C  is  spec-fled  to  create  a  contour  plot  of  the  reflection  loss  va  angle  and  fre¬ 
quency. 

1  When  option  C  is  selected,  the  41  selected  frequency  samples  are  placed  equidistantly 
on  a  logarithmic  frequency  axis.  The  sampling  in  angle  is  the  same  as  above,  but  note 
that  both  MFOU  and  HAOO  have  been  set  to  0  in  order  to  disable  the  creation  of  the 
normal  curve  plots  produced  in  case  1 . 

3  The  parameters  for  the  angle  axis  are  given  exactly  as  for  the  curve  plots  in  case  1. 
However,  for  the  logarithmic  frequency  axis  the  total  length  of  the  axis  is  not  specified. 
Instead  the  length  of  one  <  ctave  is  given  in  cm,  here  3  cm.  The  tick  marks  are  controlled 
by  the  last  parameter  in  the  line,  in  this  case  one  tick  mark  per  octave  is  plotted. 

4  The  last  line  in  the  data  file  specifies  that  contour  lines  are  to  be  plotted  for  each  1  dB 
in  the  interval  l  to  10  dB. 
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SD=  .0  m  R 


CONFR.P-P 


Grazing  angle  (degrees) 
SAFARI-FIPR  case  2. 


SD=  0.0m  RD=  0.0m 


CONFR.P-P 


V' 


U 

a 

v 

O’ 

1 40 


BH  ABOVE  10.0 

Ih|  b.0  -  io.o 

■■  1(0  -  6.0 

[  |  BELOV  2.0 


0.0  30.0  80.0  90.0 

Grazing  angle  (degrees) 

SAFARI-FIPR  case  2. 

Fig.  7.  Results  for  SAFARI-FIPR  case  2:  (a)  contour  plot  of  reflection  loss  ns 
function  of  grasing  angle  and  frequency  produced  by  the  DISSPLA  package; 
(b)  contour  plot  produced  by  UN  IRAS. 
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■  6.3.3.  SAFARI-FIPR  case  3:  Reflection  loss  at  rough  ice  layer 
The  last  example  of  using  SAFARI-FIPR  concerns  the  determination  of  the  reflection 
loss  appearing  at  a  randomly  rough  ice  layer.  This  problem  from  arctic  acoustics 
is  included  for  two  purposes.  First  of  all  it  illustrates  how  to  include  the  effect  of 
interface  roughness  in  the  calculations.  Secondly  it  is  an  example  of  a  problem  where 
the  environment  has  to  be  turned  upside  down  as  SAFARI-FIPR  always  places  the 
source  in  the  upper  halfspace. 

It  is  assumed  that  an  ice  layer  of  average  thickness  3.9  m  covers  a  water  column  with 
a  constant  sound  speed  of  1500  m/s.  The  free  surface  of  the  ice  layer  is  stochastically 
rough  with  an  rms  elevation  of  0.6  m.  The  underside  of  the  ice  has  an  rms  roughness 


Table  7 

SAFARI-FIPR  case  3a  data  file 


Data  file 

Description 

Notes 

SAFARI-FIPR  case  3a 

title 

N 

option  N  -  normal  reflection  loss 

1 

3 

number  of  layers 

0  1600  00010 

layer  1  (upper  halfspace,  water) 

O 

dl 

0  3000  1600  1.0  2.6  0.9  1.9 

layer  2  (ice  layer) 

3 

3.9  0  0  0  0  0  0.6 

layer  3  (vacuum) 

4 

1  60  60  0 

frequency  interval  1  to  50  Hz,  50  samples 

5 

10  10  1  1 

angle  interval  10°,  1  sample. 

5 

0  60  20  10 

frequency  axis  for  plot 

6 

0  0.2  12  0.1 

loss  axis  for  plot 

6 

1  Option  N  means  that  the  default  reflection  loss  is  to  be  calculated.  This  op. ion  need 
not  to  be  specified,  but  the  option  field  should  always  be  represented  by  one  line,  which 
in  this  case  could  be  empty. 

1  The  environment  has  been  turned  upside  down.  Therefore  the  upper  halfspace  is  water. 
Note  that  the  depth  and  roughness  parameters  are  dummy  for  the  upper  halfispace, 
the  depth  of  the  first  interface  and  its  roughness  data  should  be  given  together  with 
the  material  properties  of  the  first  real  layer,  i.e.  layer  2. 

3  The  roughness  of  the  water/ice  interface  is  specified  to  be  1.9  m.  Since  this  is  a 
positive  number,  the  correlation  length  is  assumed  to  be  infinite  and  has  therefore 
been  omitted.  The  scattering  theory  based  on  the  Kirchhoff  approximation  will  be 
applied  in  this  case. 

4  The  roughness  of  the  free  ice  surface  is  specified  to  be  0.6  m,  again  the  correlation 
length  is  infinite. 

8  The  frequency  interval  is  chosen  to  be  1  to  50  Hz,  and  since  the  contour  option  C  was 
not  selected,  the  50  samples  will  be  placed  equidistantly  on  a  linear  frequency  axis. 
The  calculation  is  to  be  performed  at  only  one  grazing  angle,  NANG  =  1  and  NAOU  =  1. 

a  In  contrast  to  the  case  for  the  contour  plot,  the  frequency  axis  is  here  linear,  and  the 
parameters  are  given  in  the  usual  way. 
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of  1.9  m.  Initially  the  correlation  length  is  assumed  to  be  infinite.  The  ice  has  a 
compressional  velocity  of  cc  =  3000  m/'s,  a  shear  velocity  of  c,  =  1600  m/s  and  a 
density  of  0.9  g/cm3.  The  compressional  and  shear  attenuations  are  1.0  dB/A  and 
2.5  dB/A,  respectively.  The  data  file  given  in  Table  7  is  set  up  to  calculate  the 
reflection  loss  as  function  of  frequency  at  a  fixed  grazing  angle  of  10°. 

The  resulting  plot  of  the  reflection  loss  as  function  of  frequency  is  shown  in  Fig.  8a. 
Note  the  expected  increase  in  the  loss  with  increasing  frequency.  Now  assume  that 
the  correlation  length  for  the  roughness  of  the  water /ice  interface  is  40  m,  and  not 
infinite  as  above.  Then  the  data  file  would  look  as  shown  in  Table  8.  As  can  be 
observed,  the  data  file  is  identical  to  the  one  above  except  for  the  line  stating  the 
environmental  data  for  the  ice  layer. 


Table  8 

SAFARI-FIPR  case  3b  data  file 


Data  file 


Description 


Notes 


JiFi&I-FIPK  casa  3b 

N 

3 

o  1600  00010 
0  3000  1000  1.0  2.6  0.9  -1.9  40 
3.9  0  0  0  0  0  0.8 
1  60  SO  0 
10  10  1  1 
0  60  20  10 
0  0.2  12  0.1 


title 

option  N  -  normal  reflection  loss 

number  of  layers 

layer  1  (upper  halfspace,  water) 

layer  2  (ice  layer) 

layer  3  (vacuum) 

frequency  interval  1  to  50  Hz,  50  samples 
angle  interval  10°,  1  sample 
frequency  axis  for  plot 
loss  axis  for  plot 


1 


1  By  specifying  the  negative  rms  value  of  the  roughness,  the  non-Kirchhoff  scattering 
theory  will  be  applied,  and  in  this  case  the  correlation  length  (40  m)  has  to  be  given 
as  a  last  number  in  the  same  line. 


The  resulting  reflection-loss  curve  is  shown  in  Fig.  8b.  There  is  a  significant  effect  of 
the  finite  correlation  length  in  the  actual  frequency  interval,  but  at  higher  frequencies 
the  non-Kirchhoff  result  will  approach  the  Kirchhoff  result.  More  examples  treating 
interface  roughness  can  be  found  in  [24,25]. 
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(a)  REFLECTION  COEFFICIENT 


(b)  REFLECTION  COEFFICIENT 


Fig.  8.  Results  for  SAFARl-FIPR  case  3:  (a)  reflec¬ 
tion  loss  at  rough  ice  layer  as  function  of  frequency 
at  grasing  angle  10*,  correlation  length  L  =  oo; 
(b)  reflection  loss  at  rough  ice  layer  as  function  of 
frequency  at  grasing  angle  10*,  correlation  length 
L  =  40  m. 
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7=  Running  SAFARI-FIP 

The  SAFARI-FIP  module  performs  single- frequency  calculation  of  the  total  wavefield 
at  any  number  of  depths  and  ranges  for  either  a  single  compressional  source  or  a 
vertical  phased  source  array.  In  both  cases  the  sources  may  be  either  point  sources 
or  line  sources,  yielding  axisymmetric  and  plane  fields,  respectively. 

In  the  case  of  a  single  source  in  a  fluid  medium,  the  source  strength  5W  is  normalized 
to  yield  a  pressure  of  1  Pa  at  1  m  distance  from  the  source.  For  both  a  compres¬ 
sional  source  and  a  point  force  in  a  solid  medium  the  source  strength  is  normalized 
to  yield  a  normal  stress  <rxt  of  1  Pa,  1  m  below  the  source.  The  graphical  results 
are  always  given  in  dB  and  for  the  pressure  or  normal  stress  they  therefore  directly 
correspond  to  the  standard  definition  of  the  transmission  loss.  For  the  phased  ar¬ 
rays,  the  source  strength  of  each  individual  source  is  divided  by  the  total  number 
of  sources  in  order  to  get  a  maximum  level  in  the  generated  beams  which  is  inde¬ 
pendent  of  the  number  of  sources.  It  is  common  in  both  underwater  acoustics  and 
seismology  to  measure  particle  velocities  by  means  of  geophones,  and  these  param¬ 
eters  are  therefore  included  as  optional  outputs.  Since  the  particle  velocities  are 
vector  parameters  depending  not  only  on  the  receiver  position  but  also  on  the  di¬ 
rection,  there  is  no  logical  ‘transmission  loss’  definition  as  in  the  case  of  the  scalar 
pressure.  In  the  graphical  output  these  parameters  are  therefore  given  directly  in 
dB  relative  to  1  m/s  for  a  source  strength  of  1  Pa  at  1  m  distance.  Thus  the  dB 
values  for  pressure  and  particle  velocity  will  in  general  be  offset  with  respect  to  each 
other.  For  a  plane  wave  in  a  fluid  with  density  p  and  sound  speed  cc,  the  offset  is 
determined  by  the  relation  between  pressure  p  and  particle  velocity  t>,  |p|  =  pcc\v\, 
i.e.  approximately  123  dB  for  water. 

For  all  parameters  the  output  is  available  both  as  standard  curve  plots  vs  depth  or 
range  and  as  contour  plots  showing  the  field  variation  in  both  depth  and  range. 

All  inputs  to  SAFARI-FIP  are  read  from  the  file  currently  assigned  to  the  logical  file 
FOROOl.  Before  running  the  program  the  user  has  to  assign  the  file  containing  the 
input  data  to  this  logical  name.  The  most  convenient  approach  is  to  include  this 
assignment  in  a  general  command  file  which  also  assigns  file  names  to  the  logical 
names  of  the  output  files  as  described  in  Subsect.  5.6. 

A  successful  use  of  SAFARI-FIP  is  highly  dependent  on  how  the  parameters  are 
specified,  in  particular  those  related  to  the  truncation  and  discretization  of  the 
horizontal  wavenumber  interval.  We  will  therefore  in  the  following  first  describe  the 
preparation  of  the  input  files  in  detail  and  then  outline  some  important  numerical 
considerations. 
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Tabie  9 

Parameters  of  SAFARI-FIP  input  files:  Calculation  and  environmental  parameters 


Block 

Parameter 

Units 

Limits 

I 

TITLE  : 

title  of  run 

- 

<  80  char. 

II 

optl  opt2  ...: 

output  options 

- 

<  40  char. 

III 

FREQ: 

source  frequency 

Hs 

>  0 

COFF: 

integration  contour  offset 

dB/A 

COFF  >  0 

IV 

NL: 

number  of  layers,  incl.  halfspaces 

- 

NL  >  2 

D: 

depth  of  interface 

m 

- 

CC: 

compressional  speed 

m/s 

CC  >  0 

CS: 

shear  speed 

m/s 

- 

AC: 

compressional  attenuation 

dB/A 

AC  >  0 

AS: 

shear  attenuation 

dB/A 

AS  >  0 

RQ: 

density 

g/cm3 

RO  >  0 

RG: 

rms  value  of  interface  roughness 

m 

- 

CL: 

correlation  length  of  roughness 

in 

CL  >  0 

s  r 

V 

SB: 

source  depth  (mean  for  array} 

III 

- 

NS: 

number  of  sources  in  array 

- 

NS  >  0 

DS: 

vertical  scurce  spacing 

m 

DS  >  0 

AN: 

grasing  angle  of  beam 

deg 

- 

Ik: 

array  type 

- 

1  <  IA  <  5 

FD: 

focal  depth  of  beam 

m 

FD  ^  SD 

VI 

RD1: 

depth  of  first  receiver 

m 

- 

RD2: 

depth  of  last  receiver 

m 

RD2  >  RD1 

NR: 

number  of  receivers 

- 

NR  >  0 

IR: 

plot  output  increment 

- 

IR  >  0 

VII 

CHIN: 

minimum  phase  velocity 

m/s 

CNIN  >  0 

CNAZ: 

maximum  phase  velocity 

m/s 

VIII 

NV: 

number  of  wavenumber  samples 

— 

NV  =  2m 

IC1: 

first  sampling  point 

- 

IC1  >  1 

IC2: 

last  sampling  point 

— 

IC2  <  NV 

7.1.  INPUT  FILES  FOR  SAFARI-FIP 

The  input  data  are  structured  in  13  blocks.  The  first  eight,  shown  in  Table  9, 
specify  the  title,  options,  frequency,  environmental  parameters,  source  and  receiver 
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geometry,  and  finally  the  wavenumber  sampling  parameters.  The  last  5  blocks, 
outlined  in  Table  10,  contain  axis  specifications  for  the  graphical  output.  Some  of 
these  blocks  should  always  be  included  and  others  only  if  certain  options  have  been 
specified.  The  single  blocks  and  parameters  are  described  in  detail  in  the  following. 


Table  10 

Parameters  of  SAFARI-FIP  input  files:  Plot  parameters 


Block 

Parameter 

Units 

Limits 

IX 

RMIN: 

minimum  range  on  plots 

km 

- 

MUX: 

maximum  range  on  plots 

km 

- 

ELEN: 

length  of  x-axis  for  all  plots 

cm 

ELEN  >  0 

EINC: 

distance  between  tick  marks 

km 

EINC  >  0 

X‘ 

THIN: 

minimum  transmission  loss 

dB 

- 

THAI: 

maximum  transmission  loss 

dB 

- 

TLEN: 

length  of  vertical  TL  axes 

cm 

TLEN  >  0 

TINC: 

distance  between  tick  marks 

dB 

TINC  >  0 

XI1 

DCCP: 

minimum  depth  for  plots 

m 

— 

DCLO: 

maximum  depth  for  plots 

m 

- 

DCLN: 

length  of  depth  axis 

cm 

DCLN  >  0 

DCIN: 

distance  between  tick  marks 

m 

DCIN  >  0 

XII3 

ZNIN: 

minimum  contour  level 

dB 

- 

ZNAX: 

maximum  contour  level 

dB 

- 

ZINC: 

contour  level  increment 

dB 

ZINC  >  0 

XIII4 

VLEF : 

wave  speed  at  left  border 

m/s 

— 

TEIG: 

wave  speed  at  r  ght  border 

m/s 

- 

TLEN: 

length  of  wave  t  peed  axis 

cm 

TLEN  >  0 

TINC: 

wave  speed  tick  mark  distance 

m/s 

TINC  >  0 

DTOP: 

depth  at  upper  border 

m 

- 

DTLO: 

depth  at  lower  border 

m 

- 

DTLN: 

length  of  depth  axis 

cm 

DTLN  >  0 

DTIN: 

depth-axis  tick  mark  interval 

m 

DTIN  >  0 

1  Only  for  options  i,  D,  I  and  T.  3  Only  for  options  C  and  D.  3  Only  for  option  C. 
4  Only  for  option  Z, 
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Biom  i  - 

tms 

lines: 

1 

Required: 

Always 

Syntax: 

TITLE 

Format: 

A80 

Description: 

TITLE:  Title  of  run.  Maxim  *  80  characters. 
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-  OPTIONS 


1 

Always 

optl  opt 2  opt3  . . . 

Free.  Maximum  40  alphanumeric  characters. 

Scries  of  alphanumeric  characters  specifying  the  desired 

options.  The  following  are  implemented  at  present: 

A:  Depth-averaged  transmission  loss  plotted  for  each 
of  the  selected  field  parameters  (M,  V,  H).  The 
averaging  is  performed  over  the  specified  number 
of  receivers  (Block  VI). 

C*  Range/ depth  contour  plot  for  transmission  loss. 

A  separate  contour  plot  will  be  produced  for 
every  selected  parameter  (N,  V,  H). 

D:  Plots  of  transmission  loss  vs  depth  produced 
for  the  ranges  determined  by  the  range  axis 
parameters  RMIN,  UMAX  and  RINC  (Block  IX). 

B :  Calculation  of  horisontal  particle  velocity. 

I :  Hankel  transform  integrands  are  plotted  for  each 
of  the  selected  field  parameters  («,  V,  B). 

J .  Complex  integration  contour.  The  contour 
is  shifted  into  the  upper  halfplane  by  an 
offset  controlled  by  the  input  parameter  CQPF 
(Block  ID). 

L:  Linear  vertical  source  array. 

(I ;  Calculation  of  normal  stress  <rMt  (equal  to 
negative  pressure  in  fluids). 

P:  Plane  geometry.  The  sources  will  be  line  sources 
y^jjjuteftd  of  point  sources  as  used  in  the  default 
:  cylindrical  geometry. 

T;  Transmission  Iocs  plotted  as  function  of  range  for 
each  of  the  selected  field  parameters  (N,  V,  B). 

V :  Calculation  of  vertical  particle  velocity. 

X :  If  sources  are  present  in  a  solid  medium,  these 
will  be  considered  vertical  point  forces  (line  forces 
in  plane  geometry)  rather  than  omnidirectional 
sources. 

Z :  Plot  of  velocity  profile. 
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BLOCK  m  - 

Lines: 
Required: 
Syntax: 
Format: 
Description : 


FREQUENCY 

1 

Always 
FREQ,  COFF 
Free 

FREQ:  Frequency  /  in  Hz,  at  which  calculation  has  to  be 
performed. 

COFPj  Integration  contour  offset.  Specified  in  dB/A, 
where  A  is  the  wavelength  at  the  source  depth 
$D.  Since  only  the  horizontal  part  of  the 
integration  contour  is  considered,  this  parameter 
should  not  be  chosen  so  large  that  the  amplitudes 
at  the  ends  of  the  integration  interval  become 
significant.  In  lossless  cases  too  small  values  will 
give  sampling  problems  at  the  normal  modes  and 
other  singularities.  For  intermediate  values*  the 
result  is  independent  of  the  choice  of  COFF,  but  a 
good  value  to  choose  is  one  that  gives  60  dB  at 
the  longest  range  considered  in  the  FFT,  i.e. 

COFF  =  60cc(z4)//iw 
where  the  maximum  FFT  range  is 

rm»*  =  N"W/ /(cro|n  —  Cn,1,,). 

This  value  is  the  default  which  is  applied  if  COFF 
is  specified  as  0.0. 
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ENVIRONMENTAL  DATA 
>  3 

Always 


D<1)  CCCl)  CS(1)  AC(1)  AS(1)  ROil)  RG(1)  CL(1) 

D(2)  CC(2)  CS(2)  AC(2)  AS(2)  R0<2)  RG(2)  CI.(2> 

• 

D(NL)CC(NL)CS(NL)AC(NL)AS(NL)RO(NL)RG(NL)CL(HL) 


Number  of  layer*,  including  the  upper  and  lower 
halfspaces.  These  should  always  be  included,  even 
in  cases  where  they  are  vacuum. 

Depth  z  in  m  of  upper  boundary  of  layer  or 
halfspace.  The  reference  depth  can  be  choosen 
arbitrarily,  and  DO  is  allowed  to  be  negative. 

For  layer  no.  1,  i.e.  the  upper  halfspace,  this 
parameter  is  dummy. 

Velocity  ce  of  compressional  waves  in  m/s. 

If  specified  as  0.0,  the  layer  or  halfopace  is  a 


CC(): 


CSC) : 


vacuum. 


ACO: 


AS<): 


Velocity  c,  of  shear  waves  in  m/s.  In  order  to  be 
physically  meaningful,  it  is  required  that  c,  < 
V0.75cc.  If  specified  as  0.0,  the  layer  or  hedfspace 
is  fluid.  If  c,  <  0,  it  represents  the  compressional 
velocity  at  the  bottom  of  the  actual  layer,  which 
is  treated  as  a  fluid  with  i/c(.*)*  varying  linearly 
with  depth. 

Attenuation  jc  of  compressional  waves  in  dB/A. 
If  the  layer  is  fluid,  and  ACO  is  specified  as  0.0, 
then  an  empirical  water  attenuation  is  used. 
Attenuation  7,  of  shear  waves  in  dB/A.  For  the 
attenuations  to  be  physically  meaningful,  it  is 
required  that 


*5  ft)' 
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BLOCK  IV  -  (Continued) 

Description: 

R0():  Density  p  in  g/cm*. 

RG()  :  RMS  roughness  of  interface  in  m.  RG(1)  is 

dummy.  If  RG  >  0  the  KirchhofT  approximation 
is  used  for  the  scattering  loss,  i.e.  CL  -♦  oo. 

This  is  computationally  faster  than  doing  the 
full  scattering  theory,  but  it  is  only  applicable  for 
small  roughnesses  and  large  correlation  lengths. 
Tb  invoke  the  full  non-Kirchhoff  scattering  RG 
must  be  specified  equal  to  the  negative  rms 
roughness,  i.e.  RG  <  0. 

CL() :  Correlation  length  of  roughness  in  m.  This 
parameter  is  only  required  if  non-Kirchhoff 
scattering  has  been  selected,  i.e.  RG  <  0. 
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BLOCK  V  ~ 

SOURCE  DATA 

Lines: 

1 

Required: 

Always 

Syntax: 

Option  L:  SD  NS  OS  AN  IA  FO 

format; 

else:  SD 

Free 

Description: 

SD :  Source  depth  x,  in  m.  If  option  L  has  been 

specified,  SD  defines  the  mid-point  of  the  vertical 
source  array 

NS :  Number  of  sources  in  the  array. 

DS:  Source  spacing  in  m. 

AN :  Specifies  the  nominal  grazing  angle  of  the 

generated  beam  in  degrees.  AN  >  0  corresponds 
to  downward  propagation. 

IA  s  Array  type 

(1]  Rectangular  weighted  array 

[2]  Hanning  weighted  array 

(3]  Hanning  weighted  focused  array 

[4]  Gaussian  weighted  array 

[5]  Gaussian  weighted  foci,  sed  array 

FD:  Focal  depth  in  m  for  an  array  of  type  3  or  5. 
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B^OCK  VI  -  RECEIVER  DATA 


Lines: 

Required: 

Syntax: 

Format: 

Description: 


1 

Always 

RD1  RD2  Nk  IR 

Pres 

RD1 :  Receiver  depth  in  m.  If  the  calculations  are  to  be 
performed  for  mote  them  one  receiver,  i.e.  NR  >  1  , 
then  RD1  specifies  the  uppermost  receiver. 

1U»2 :  If  NR  >  1  this  parameter  determines  the  depth 
of  the  lowermost  receiver.  If  NR  =  1,  then  this 
parameter  is  dummy. 

NR:  N’lmber  of  receivers  in  depth. 

IR:  If  options  T  or  I  have  been  specified,  plots  will 

be  generated  for  each  IR  receiver  depth.  This 
parameter  is  useful  if  option  C  or  A  have  been 
choc  sen,  and  integrands  or  transmission  loss 
curves  are  wanted  only  for  a  limited  number  of 
depths. 
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BLOCK  VII  - 

PHASE  VELOCITIES 

liner: 

1 

Required: 

Always 

Syntax: 

CHIN  CMAX 

Format: 

Description: 

Free 

Minimum  phase  velocity  cW|n  in  m/s,  Determines 

CHIM: 

the  upper  limit  of  the  truncated  horizontal 
wavenumber  space: 

=  f l  Cmin- 

CMAX: 

Maximum  phase  velocity  cmax  in  m/s. 

Determines  the  lower  limit  of  the  truncated  ! 
horizontal  wavenumber  space: 

^min  =  f/Cfntx- 

In  plane  geometry  (option  P)  CMAX  may  be 
specified  as  negative.  In  ttus  case,  the  negative 
wavenumber  spectrum  will  be  included  with 
*mi«  =  “*W*>  yielding  correct  solution  also 
at  zero  range.  Contour  offset  (option  J  )  is  not 
allowed  in  this  case. 
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BLOCK  vm  - 

WAVENUMBER  SAMPLING 

Lines: 

1 

.Required; 

Always 

Syntax: 

NW  IC1  IC2 

...  Bomtkt: 

Free 

Description: 

HW: 

Number  of  sampling  points  in  wavenumber  space; 
Should  be  an  integer  power  of  2,  i.e.  NW  =  2**. 

The  sampling  points  are  placed  equidistantly  in 
the  truncated  wavenumber  space  determined  by 
CHIN  and  CMAX 

ICt: 

Number  of  the  first  sampling  point,  where  the 
calculation  is  to  be  performed.  If  IC1  >  1, 
then  the  Hankel  transform  is  seroed  for  sampling 
points  1,2... IC1  -  1,  and  the  discontinuity  is 
smoothed. 

IC2t 

Number  of  the  last  sampling  point  where  the 
calculation  is  to  be  performed.  If  1C2  <  NW, 
then  the  Hankel  transform  is  zeroed  for  sampling 
points  ICS  +  1 , . . ,  NW,  and  the  discontinuity  is 
smoothed  by  Hermite  polynomial  extrapolation. 

I 

\ 
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BLOCK  JX  ~ 

PLOT  DAM  /  RANGE  AXES 

Lines: 

1 

Required; 

Always 

Syntax : 

RUIN  RHAX  RLEN  RING 

Pbrm&t: 

Eree 

Description: 

RMIN :  Minimum  range  rn,^  in  km.  Determines  the 

starting  value  for  the  x-axis  of  transmission  loss 

and  contour  plots.  RMIN  is  also  the  first  range  for 
which  transmission  loss  vs  depth  will  be  plotted  if 
option  D  is  selected. 

RMAI:  Maximum  range  rm%x  in  km.  This  parameter 

only  defines  the  end  of  the  x-axis  of  transmission 
loss  and  contour  plots.  The  range  interval  of 

calculation  is  defined  by  rmjn  and  the  truncation 

and  discretisation  of  the  wavenumber  space. 

RLEN :  Length  of  x-axis  in  cm  of  all  plots,  i.e.  also  the 
length  of  the  transmission  loss  axis  if  option  0  is 

selected. 

RINCt  Range  interval  between  tick  marks  on  the  x- 
axes  of  transmission  loss  and  contour  plots. 

Also  determines  the  range  interval  between 
plots  of  transmission  loes  vs  depth  for  option  0. 
Fbr  integrand  plots  all  scaling  and  labelling  is 

automatic. 
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BLOCK  X  ~  PLOT  DATA  /  TRANSMISSION  LOSS  AXES 

Line*:  1  for  each  field  parameter  (M,  ¥,  B). 

Required:  If  option*  A,  D,  I  or  T  were  specified. 

Syntax:  TRIM  TMAX  TLEN  TIBC 

jFhrmat:  IVee 

Description: 

THIN:  Minimum  transmission  loss  in  dB.  For  plots 
of  transmission  loss  vs  range  (options  A  and 
T),  THIN  determines  the  upper  border.  In  the 
case  of  transmission  loss  vs  depth  (option  0),  it 
determines  the  right  border  of  the  plot. 

TMAX :  Maximum  transmission  loss  in  dB.  Determines 
the  lower  border  (for  options  A  and  T)  or  left 
border  (for  option  0)  of  transmission  loss  plots. 

It  should  be  noted,  that  the  reference  for  the 
particle  velocities  is  1  m/s,  when  the  reference 
pressure  is  1  Pa  (1  N/m2)  at  1  m  distance 
from  the  source.  If  option  L  has  been  specified, 
i.e.  more  than  one  source  is  present,  then  the 
reference  pressure  of  each  source  is  divided  by 
the  number  of  sources.  This  ensures  that  the 
order  of  magnitude  of  the  field  parameters  are 
independent  of  the  number  of  sources  in  the 
array. 

TLBH :  Length  of  y-axes  in  cm  for  transmission  loss  and 
integrand  plots.  For  option  0,  the  transmission 
loss  is  represented  by  the  axis  and  the  length  is 
;;.v':;::.::th«re*bre  determined  by  RLEN  . 

TIICs  Number  of  dB  between  tick  marks  on 

transmission  loss  axes.  Scaling  and  labeling  of 
integrand  plots  is  automatic. 
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BLOCK  XI 

PLOT  DATA  /  DEPTH  AXES 

Him: 

1 

Reqtured: 

If  options  C  or  D  were  specified. 

Syntax; 

DCUP  DCLO  DCIN  DCIN 

format: 

■  FVee 

Description: 

OCOP:  Minimum  depth  in  m  for  both  contour  plots  and 

plots  of  transmission  loss  vs  depth. 

DCLO :  Maximum  depth  in  m  for  both  contour  plots  and 

plots  of  transmission  loss  vs  depth. 

DCLN:  Length  in  cm  of  depth  axes, 

DCIN:  Difference  in  m  between  tick  marks  on  depth 

U. _ _ _ - _ i. - 

axes. 

BLOCK  xir  ~ 

PLOT  DATA  /  CONTOUR  LEVELS 

Lines: 

1  for  each  selected  parameter  (N,  V,  H). 

Required:  .  ' 

If  option  C  was  specified. 

Syntax: 

ZMIN  ZMAZ  ZINC 

'■■format:-  ,■ 

Btee 

'  Description:  -I 

ZMIN :  Lower  limit  in  dB  of  transmission  loss  interval  for 
which  contour  lint*  are  to  be  plotted. 

ZHAX  i  Upper  limit  in  dB  of  transmission  loss  interval  for 
which  contour  lines  are  to  be  plotted. 

ZINC;  Difference  in  dB  between  contours. 
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BLOCK  xm  ~ 

BLOT  DATA  /  VELOCITY  PROFILE 

Ltow 

2 

Required: 

If  option  Z was  specified  in  Block  II. 

Synt&x: 

VUBF  VRIG  VLEN  VIMC 

DVOP  DVLO  DVLN  DVIN 

FbmtA: 

Pnte 

UemCtipuOa: 

VLEF :  Left  border  value  of  velocity  in  m/s. 

VRIG:  Right  border  value  of  velocity  in  m/s. 

VUKN:  Length  of  velocity  axis  in  cm. 

VIRC:  Distance  between  tick  marks  in  m/s. 

DVUP;  Depth  at  upper  border  of  plot  in  m. 

DVLO :  Depth  at  lower  border  of  plot  in  m. 

DVLN :  Length  of  depth  axis  in  cm. 

DVIN:  Distance  in  m  between  tick  marks  on  depth  axis. 
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7.2.  NUMERICAL  CONSIDERATIONS 

When  SAFARI-FIP  is  to  be  applied  to  a  propagation  problem,  the  first  difficulty  that 
arises  is  related  to  the  discretization  of  the  environmental  model.  When  modelling 
real  data  collected  in  an  area  where  only  a  few  environmental  data  are  known, 
it  is  generally  necessary  to  re-run  the  code  several  times  in  order  to  lit  a  discrete 
stratified  model;  on  the  other  hand  the  accuracy  requirements  will  usually  be  limited, 
and  there  will  be  no  purely  numerical  considerations  involved.  However,  the  choice 
of  environmental  model  is  important  in  cases  where  the  SAFARI  results  have  to  be 
compared  to  those  obtained  by  other  numerical  models,  which  is  very  common  in 
relation  to  the  verification  of  new  algorithms.  Many  codes,  like  the  normal  mode  and 
parabolic  equation  models,  assume  a  linearly  varying  sound  velocity  profile  between 
the  sampling  points,  whereas  SAFARI  assumes  either  isovelocity  layers  or  a  sound 
speed  varying  non-line&rly  between  the  sampling  points,  Eq.  (5).  It  is  therefore 
necessary  for  the  number  of  profile  sampling  points  to  be  so  large  that  the  difference 
between  the  interpolation  techniques  has  no  significant  effect  on  the  calculated  sound 
field.  It  is  not  possible  to  give  any  specific  rules  in  this  regard  since  the  effect  depends 
on  many  factors  such  as  frequency,  water  depth  and  the  receiver  range  of  interest. 
But  if  isovelocity  layers  are  used  to  represent  the  profile,  a  layer  thickness  of  J-A  will 
usually  be  sufficient.  For  the  inhomogeneous  layers,  however,  the  thickness  can  be 
chosen  to  be  much  larger.  In  any  event  the  profile  discretization  should  be  controlled 
by  checking  the  convergence. 

The  actual  choice  of  parameters  controlling  the  wavenumber  integration  is  far  more 
critical.  There  are  essentially  5  parameters:  CHIN,  CMAX,  NW,  IC1  and  IC2,  but 
they  cannot  be  chosen  independently  from  other  parameters,  in  particular  the  range 
interval  of  interest  (Block  IX),  because  of  the  wrap-around  problem.  Further,  both 
source  and  receiver  depths  influence  the  choice  of  integration  parameters,  e.g.  in 
the  waveguide  problem  the  presence  of  the  normal  modes  will  require  much  higher 
sampling  rate  if  the  source  and  receiver  are  inside  the  waveguide  than  if  one  of  them 
is  buried  in  the  bottom.  There  are  no  general  rules  which  can  be  used  to  automate 
the  wavenumber  integration.  The  only  way  to  obtain  accurate  results  is  to  change 
the  integration  parameters  until  convergence  is  obtained. 

Clearly,  when  mak'ng  parameter  studies,  convergence  tests  sure  not  required  for  every 
small  change  in  the  environmental  parameters.  However,  convergence  tests  must 
be  done  for  at  least  one  characteristic  example  before  proceeding  with  a  complete 
propagation  study.  A  user  with  a  reasonable  knowledge  of  the  physics  of  waveguide 
problems  will  -  after  gaining  some  experience  -  be  able  to  determine  the  proper 
parameters  relatively  quickly.  All  problems  are,  however,  not  equally  easy  to  predict, 
and  it  is  therefore  advisable  -  also  for  the  experienced  user  -  to  use  the  following 
guidelines  for  estimating  the  integration  parameters  for  every  new  application  of 
SAFARI-FIP. 

1.  Select  the  horizontal  phase  velocities  CHIN  and  CMAX  such  that  all  signifi¬ 
cant  wave  phenomena  are  included,  i.e.  CMIN  should  be  chosen  to  be  10-20% 
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smaller  than  the  smallest  wave  speed  in  the  problem  and  CMAX  should  be  a 
very  large  number,  e.g.  CMAX  =  1E8.  If  evanescent  waves  are  known  to  be 
present  CHIN  may  have  to  be  even  smaller. 

2.  Select  a  few  characteristic  receiver  depths. 

3.  Select  a  reasonably  small  number  of  sampling  points  NW,  and  run  the  code 
with  option  1  to  obtain  plots  of  the  integration  kernels. 

4.  Repeat  steps  1  to  3  with  decreasing  CHIN  until  all  significant  wavenumber 
components  are  within  the  integration  interval.  Then  read  the  wavenumber 
limits  Jbmin  and  Jfem«  from  the  plot  or  calculate  them  directly  as  fcTO|n  = 
2ir//cmM  and  =  2x//cm in- 

5.  Select  the  range  interval  (RM1N,  RMAX)  of  interest,  and  convert  from  km  to  m 

to  obtain  rmi„  and  rm„.  For  the  FFT  integration  to  ‘reach’  rmM  the  nec¬ 
essary  number  of  sampling  points  must  be  larger  than  Nmin  =  - 

fcmin)/2ir}.  Choose  the  smallest  M  for  which  NW  =  2M  >  Nmia.  Set  IC1  =  1 
and  IC2  =  NW. 

6.  If  the  integration  kernels  plotted  out  in  step  3  are  very  ‘peaky’  due  to  low  loss 
in  the  waveguide,  select  option  J  and  specify  COFF  =  0  to  invoke  the  default 
contour  offset.  The  experienced  user  can  choose  a  specific  value  of  COFF  as 
described  in  Block  III  above. 

7.  Now  select  option  T  and  calculate  the  transmission  loss. 

8.  Double  the  number  of  sampling  points  NW  (remember  to  double  IC2  as  well) 
and  re-calculate. 

9.  Repeat  step  8  until  a  stable  transmission  loss  curve  is  obtained. 

10.  Now  the  number  of  receivers  can  be  increased  and  the  other  options  C,  A 
and  Z  can  be  specified,  as  desired.  Remember  to  change  the  plot  parameters 
accordingly,  see  Table  10. 

This  is  the  standard  procedure  which  is  only  applicable  if  the  number  of  sampling 
points  NW  is  reasonably  small,  at  least  smaller  than  the  m;  isim  allowed  for  the 
actual  installation.  If  this  is  not  the  case,  there  are  two  different  ways  to  proceed.  If 
the  computation  time  is  not  important  find  the  actual  computer  allows  it,  a  larger 
version  can  be  installed,  and  the  procedure  above  is  repeated  or  continued.  In  many 
cases,  however,  the  computation  time  will  become  unacceptable  using  this  approach, 
and  the  second  possible  approach  has  to  be  taken.  This,  however,  requires  a  little 
more  knowledge  concerning  the  physics  of  the  actual  problem. 

As  discussed  earlier,  each  horizontal  wavenumber  component  k  corresponds  to  a 
specific  conical  or  plane  wave  propagating  at  grazing  angle  0,  where  Jfc  =  km  cos  0, 
km  being  the  medium  wavenumber  of  layer  m.  Therefore  the  small  wavenumbers 
correspond  to  steep  propagation,  and  k  =  km  corresponds  to  horizontal  propagation. 
It  is  well  known  for  the  waveguide  problem  that  the  waves  propagating  at  grazing 
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angles  larger  than  the  critical  angle  9C  at  the  bottom  will  suifer  a  large  attenuation 
for  every  bottom  bounce  and  will  therefore  yield  ^n  insignificant  contribution  to 
the  toted  field  at  longer  ranges.  This  so-called  continuous  part  of  the  wavenumber 
spectrum  0  <  Jfe  <  Jfec  can  therefore  be  ignored  for  long  range  propagation.  This  is  in 
fact  the  approximation  made  in  most  normal  mode  models  where  only  the  discrete 
part  Jfec  <  Jfe  <  Jfem  of  the  spectrum  is  considered.  In  SAFARI-FIP  the  wavenumber 
spectrum  can  be  reduced  by  specifying  a  smaller  value  of  CMAX.  This  will  obviously 
lead  to  a  smaller  wavenumber  spectrum,  and  thus  to  larger  ranges  covered  by  the 
FFT  integration.  Due  to  the  fact  that  the  wavenumbi  *  spectrum  is  a  continuous 
function  of  Jfe,  an  arbitrary  truncation  point  may  give  ri  to  numerical  artifacts  such 
as  WTap-around.  The  truncation  is  therefore  usually  determined  in  the  following 
way: 

1.  Carry  out  steps  1  to  2  above. 

2.  By  inspection  of  the  integration  kernels,  choose  a  truncation  point  in  the  con¬ 
tinuous  spectrum  where  the  amplitude  is  small.  Calculate  the  orresponding 
phase  velocity  ct.  If  only  a  single  receiver  depth  is  required  choose  CMAX  =  ct 
and  proceed  to  next  step.  If  more  receiver  depths  are  involved  a  universal 
truncation  point  with  vanishing  amplitude  can  usually  not  be  found.  In  these 
cases  and  when  the  kernel  has  significant  amplitude  in  the  whole  continuous 
spectrum,  choose  CMAX  somewhat  larger  than  ct.  Then  use  IC1  to  specify 
the  first  sampling  point  where  the  calculation  should  start,  and  SAFARI  will 
automatically  taper  the  discontinuity.  The  same  procedure  is  of  course  ap¬ 
plicable  if  the  wavenumber  spectrum  is  truncated  at  the  large  wavenumber 
end,  which  requires  a  proper  choice  of  CMIN  and  IC2. 

3.  Now  perform  steps  5  to  8  above. 

4.  If  determination  of  both  short  and  long  range  propagation  is  required,  cal¬ 
culate  the  full-spectrum  solution  to  a  reasonable  and  convenient  range  as 
described  above.  Then  compare  the  solutions  to  ensure  that  they  overlap.  If 
this  is  not  the  case  -  which  will  only  very  rarely  happen  --  change  the  trun¬ 
cation  point  to  include  more  of  the  continuous  spectrum  and  re-calculate  the 
long-range  field  until  there  is  an  overlapping  region  of  identical  results. 

To  the  inexperienced  user  the  procedure  outlined  here  may  seem  very  cumbersome, 
but  as  mentioned  above  the  test  of  convergence  is  not  needed  for  every  small  change 
in  the  parameters,  and  after  a  few  trials  it  will  probably  seem  a  very  logical  proce¬ 
dure. 
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7.J.  SAFARI-FIP  EXAMPLES 

Using  basically  the  same  environmental  model  as  for  the  SAFARI-FIPR  examples,  we 
will  here  demonstrate  the  different  features  of  SAFARI-FIP  by  a  series  of  examples. 
These  have  been  chosen  such  that  also  the  numerical  considerations  involved  ir  the 
wavenumber  integration  are  illustrated. 

In  all  cases  the  water  depth  will  be  100  m  and  the  ocean  bottom  is  assumed  to  have 
the  properties  given  in  Table  4. 


■  7.3.1.  SAFARI-FIP  case  1:  Seismic  interface  wave  propagation 

In  the  first  example  we  will  calculate  the  transmission  loss  at  5  Hz  for  a  source  depth 
of  95  m  and  a  receiver  depth  of  100  m  and  for  ranges  out  to  5  km.  This  is  done  by 
specifying  the  data  file  given  in  Table  1 1 . 

After  running  SAFARI-FIP  with  this  data  file,  the  plot  program  FIPPLOT  will  pro¬ 
duce  the  plots  shown  in  Fig.  9.  The  kernel  in  the  wavenumber  integral  representation 
of  the  normal  stress  at  depth  100  m  is  shown  in  Fig.  9a.  As  can  be  observed,  there 
are  two  distinct  peaks.  The  large  peak  at  a  wavenumber  of  0.07  m-1  is  the  funda¬ 
mental  interface  mode.  The  second  peak  at  a  wavenumber  of  0.017  m-1  is  a  virtual 
mode  associated  with  propagation  in  the  water  column.  Figure  9b  show*  the  calcu¬ 
lated  transmission  loss,  clearly  displaying  an  interference  pattern  produced  by  the 
interface  mode  and  the  virtual  mode.  It  should  be  pointed  out  that  the  integra¬ 
tion  contour  has  not  been  offset  in  this  case  because  no  discrete  normal  modes  are 
present.  Had  the  frequency  been  chosen  a  little  higher,  however,  the  virtual  mode 
would  have  moved  into  the  discrete  spectrum  and  given  rise  to  a  very  sharp  peak 
in  the  integration  kernel.  In  that  case  option  J  could  conveniently  be  selected,  but 
this  is  dealt  with  in  the  following  example. 

As  can  be  observed  from  Fig.  9a,  the  amplitude  of  the  integration  kernel  is  not 
insignificant  at  the  truncation  point  defined  by  the  parameter  IC2.  However,  due  to 
the  tapering  of  the  discontinuity  and  the  fact  that  the  exponential  function  is  rapidly 
varying  at  the  large  wavenumbers,  the  influence  of  the  truncation  is  insignificant. 
The  significance  of  the  truncation  can  of  course  be  checked  by  increasing  the  value 
of  IC2  and  comparing  the  results. 
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Table  11 

SAFARI-FIP  case  1  data  file 


Data  file 

Description 

Notes 

SAFiftl-FIP  case  1 

title 

NIT 

options 

1 

5.0 

frequency  5  Hs 

4 

number  of  layers 

0  0  0  0  0  0  0 

layer  1  (upper  halfspace,  vacuum) 

2 

0  1S00  00010 

layer  2  (isovelocity  water  column) 

2 

100  1600  400  0.2  0.5  1.8  0 

layer  3  (silt  sediment) 

120  1800  600  0.1  0.2  2.0  0 

layer  4  (sand  sub-bottom) 

95 

source  depth  95  m 

100  100  1  1 

receiver  depth  100  m 

100  1E8 

phase  velocities  CHIN  and  CNAX 

3 

2348  1  1000 

wavenumber  sampling  parameters 

4 

0.0  6.0  20  1.0 

range  axis  parameters 

20  80  12  10 

transmission-loss  axis  parameters 

1  Option  N  indicates  that  the  normal  stress  (=  -  pressure)  should  be  calculated.  The  I 
option  will  produce  a  plot  of  the  depth-dependent  Green’s  function,  i.e.  the  kernel  iu 
the  integral  representation  for  the  normal  stress  at  the  selected  receiver  depth.  Option 
T  will  generate  a  plot  of  the  transmission  loss  as  function  of  range. 

3  The  sea  surface  is  introduced  by  including  a  vacuum  upper  halfspace,  and  the  sea 
surface  is  chosen  as  the  origin  for  the  depth  axis. 

3  In  order  to  include  the  propagation  directions  close  to  vertical,  CHAI  has  been  set 
to  a  very  large  number  such  that  ss  0.  Since  the  source  is  close  to  the  sea¬ 
bed  and  the  frequency  is  relatively  low,  seismic  interface  waves  will  be  excited.  CHIN 
should  therefore  be  chosen  somewhat  smaller  than  the  smallest  wave  speed  in  the 
problem,  i.e.  the  sediment  shear  speed  400  m/s.  A  value  of  250  m/s  would  probably  be 
sufficiently  low.  However,  due  to  the  link  between  the  wavenumber  and  range  sampling, 
Eq.  (115),  this  would  lead  to  a  relatively  course  range  sampling,  Ar  ~  cmu  //  —  50  ra. 
To  reduce  this  to  20  m,  CHIN  is  specified  as  100  m/s. 

*  Tb«  number  of  wavenumber  sampling  points,  NW,  is  specified  as  2046  =  2n,  but  since 
no  wave  phenomena  are  expected  with  phase  velocities  smaller  than  200  m/s,  the 
depth-dependent  Green’s  function  is  only  calculated  for  the  first  1000  sampling  points 
(IC1  =  1,  IC2  =  1000),  automatically  invoking  the  discontinuity  tapering. 
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(b)  TRANSMISSION  LOSS 


Fig.  9.  Results  for  SAFARI-FIP  case  1:  (a)  depth- 
dependent  Green’s  function  as  function  of  horison- 
tal  wavenumber;  (b)  calculated  transmission  loss  as 
function  of  range. 
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■  7.3.2.  SAFARI-FIP  case  2:  Normal  mode  propagation 

As  a  second  example,  a  propagation  problem  is  chosen  where  the  interface  wave  is 
unimportant.  Because  of  the  characteristic  exponential  amplitude  decay  away  from 
the  guiding  boundary,  the  interface  wave  is  insignificant  when  either  the  source  or 
the  receiver  is  far  away  from  the  sea-bed  in  terms  of  wavelength.  Typically,  the 
contribution  to  the  total  held  by  the  interface  wave  is  insignificant  if  the  source 
is  one  wavelength  away  from  the  sea-bed  and  if  discrete  modes  are  present  in  the 
water  column.  We  will  therefore  place  a  30  Hz  source  at  mid- water  depth  and  again 
calculate  the  transmission  loss  out  to  5  km  range  for  a  receiver  at  the  ocean  bottom. 
To  demonstrate  how  a  velocity  profile  in  the  water  column  is  specified,  it  is  assumed 
that  the  sound  speed  has  been  measured  to  be  1500  m/s  at  the  surface,  1480  m/s 
at  30  m  depth,  and  1400  m/s  at  the  sea-bed.  The  data  file  is  set  up  as  shown  in 
Table  12. 

The  results  are  shown  in  Fig.  10.  Figure  10a  shows  the  integration  kernel  with  the 
velocity  profile  plot  inserted.  Two  discrete  modes  are  apparent,  and  a  third  mode  is 
just  around  cut-off  at  the  medium  wavenumber  for  compressional  waves  in  the  sub¬ 
bottom,  k  =  2sr//l800  ~  0.105.  The  resulting  transmission  loss  plotted  in  Fig.  10b 
is  clearly  dominated  by  the  interference  pattern  produced  by  the  2  discrete  modes 
for  ranges  larger  than  1  km. 
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Table  12 

SAFARI-FIP  case  2  data  file 


Data  file 

Description 

Notes 

SiFUI-FIP  case  2 

title 

N  I  T  J  Z 

options 

1 

30  0 

frequency  and  contour  offset 

2 

6 

number  of  layers 

0  0  0  0  0  0  0 

layer  1  (upper  halfspace,  vacuum) 

0  1600  -1480  0010 

layer  2  (water  layer) 

3 

30  1480  -1490  0  0  1  0 

layer  3  (water  layer) 

3 

100  1000  400  0.2  0.6  1.8  0 

layer  4  (silt  sediment) 

120  1800  600  0.1  0.2  2.0  0 

layer  5  (sand  sub-bottom) 

"0 

source  depth  50  m 

100  100  1  1 

receiver  depth  100  m 

TOO  1B8 

phase  velocities  CHIN  and  CHAI 

4 

1024  1  612 

wavenumber  sampling  parameters 

5 

0.0  6.0  20  1.0 

range  axis  parameters 

20  80  12  10 

transmission-loss  axis  parameters 

1460  1660  10  26 

velocity  axis  for  profile  plot 

1 

0  100  10  20 

depth  axis  for  profile  plot 

1 

1  The  first  3  options  are  the  same  as  specified  in  case  1,  but  here  option  J  has  been 
specified  to  invoke  the  contour  offset  since  discrete  inodes  are  present  at  this  frequency. 
Further,  option  Z  will  generate  a  plot  of  the  velocity  profile  in  the  water  column.  The 
associated  axis  specifications  are  given  at  the  end  of  the  data  file. 

1  The  frequency  is  specified  as  30  Hz  and  the  default  contour  offset  will  be  applied  since 
COFF  =  0. 

3  By  specifying  a  negative  vilue  for  the  shear  speed,  SAFARI  will  interpret  the  number 
as  the  negative  of  the  compressional  speed  at  the  bottom  of  the  layer,  i.e.  at  the  next 
interface.  The  solution  obtained  assumes  a  linear  depth  variation  of  1  /c\  within  the 
layer. 

4  Although  it  can  be  demonstrated  that  the  part  of  the  spectrum  with  phase  velocities 
smaller  than  s=s  1400  m/s  yields  no  significant  contribution  to  the  total  field,  a  relatively 
small  value  of  CHIN  has  been  selected.  This  is  again  done  in  order  to  obtain  a  reasonable 
range  sampling.  In  the  present  case  Ar  ~  cm\n/ f  =  23.3  m. 

*  As  in  the  former  case,  by  specifying  IC2  =  512  the  insignificant  second  half  of  the 
wavenumber  spectrum  is  not  calculated. 
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(b)  TRANSMISSION  LOSS 


Fig.  10.  Results  for  SAFARI-FIP  ca^e  2:  (a)  depth- 
dependent  Green’s  function  as  function  of  horison- 
tal  wavenumber  with  velocity-profile  plot  inserted; 
(b)  calculated  transmission  loss  as  function  of  range. 
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■  7.3.3.  SAFARI-FIP  case  3:  Transmission  Joss  vs  depth  and  range 

Next,  we  take  exactly  the  same  problem  considered  in  case  2,  but  calculate  the 
transmission  loss  at  several  receiver  depths,  both  in  the  water  column  and  in  the 
sediment,  and  produce  a  contour  plot  of  the  loss  as  a  function  of  range  and  depth. 
Further,  we  produce  a  plot  of  the  transmission  loss  averaged  over  depth  as  a  function 
of  range.  The  data  file  is  given  in  Table  13. 

The  resulting  plot  of  the  depth- averaged  transmission  loss  is  shown  in  Fig.  11a.  Note 
that  the  averaging  removes  the  interference  pattern  seen  at  the  individual  receiver 
depths,  Fig.  10b.  The  contour  plot  produced  by  the  CONTUR  program  is  shown  in 
Fig.  lib.  The  periodic  pattern  for  ranges  longer  than  a  couple  of  kilometres  is  typical 
for  a  2-mode  propagation  problem.  The  actual  contouring  grid  size  is  indicated  by 
the  small  ‘box’  at  the  upper  left  corner  of  the  plot. 

Figure  12  shows  two  out  of  a  total  of  5  plots  of  transmission  loss  t>s  depth  produced 
by  option  D. 
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T&ble  13 

SAFARI-FIP  cue  3  data  file 


Data  file 

Description 

Notes 

SiFAU-FIP  case  3 

title 

1  C  A  D  J 

options 

1 

30  0 

frequency  and  contour  offset 

6 

number  of  layers 

0  0  0  0  0  0  0 

layer  1  (upper  halfspace,  vacuum) 

0  1600  -1480  0  0  1  0 

layer  2  (water  layer) 

30  1480  -1490  0010 

layer  3  (water  layer) 

100  1600  400  0.2  0.5  1.8  0 

layer  4  (silt  sediment) 

120  1800  600  0.1  0.2  2.0  0 

layer  5  (sand  sub-bottom) 

50 

source  depth  50  m 

0.1  120  41  40 

receiver  depths  0.1-120  m 

2 

1350  1E8 

phase  velocities  CHIN  and  CHAX 

3 

1024  1  960 

wavenumber  sampling  parameters 

4 

0.0  6.0  20  1.0 

range  axis  parameters 

5 

20  80  12  10 

transmission-loss  axis  parameters 

0  120  12  20 

depth  axis  for  contour  plot 

40  70  6 

contour  levels  in  dB 

1  Option  C  will  create  a  contour  plot  of  the  transmission  loss  u  a  function  of  depth  and 
range,  whereu  option  k  will  calculate  the  depth-averaged  transmission  loss  over  the 
specified  number  of  receiver  depth.  Option  D  will  generate  plots  of  transmission  loss 
vm  depth  at  the  ranges  defined  by  the  range  axis  parameters,  see  note  5  below. 

1  41  receivers  will  be  placed  equidistantly  in  the  depth  interval  0.1-120  m.  Note  that  the 
first  receiver  is  not  placed  on  the  surface  where  the  field  is  known  to  vanish  and  thus 
hu  an  undefined  dB  level.  The  last  parameter  in  this  line  is  IR  =  40.  This  parameter 
is  dummy  in  the  present  cue;  but  if  option  I  or  T  had  been  specified,  an  integrand  or 
transmission  loss  plot  would  be  created  for  every  40  receiver  depth,  i.e.  here  for  the 
first  and  lut  depth  only. 

3  Compared  to  cue  2,  the  wavenumber  interval  is  smaller  here,  since  a  very  fine  range 
sampling  is  not  crucial  for  the  contour  plots  u  a  consequence  of  the  CONTUR  program 
performing  a  smooth  interpolation  between  the  data  points.  The  present  choice  of  CHIN 
and  CHAX  translates  into  a  range  step  of  Ar  ~  45  m. 

4  The  tapering  L.  again  activated  by  'pecifying  IC2  =  950. 

*  The  range  axis  parameters  are  applied  to  both  the  plot  of  the  depth-averaged  trans¬ 
mission  lou  and  the  contour  plot.  Further,  they  determine  the  ranges  for  which  trans¬ 
mission  loss  vs  depth  will  be  plotted  (option  D).  Thus,  these  plots  will  be  produced  at 
the  ranges  corresponding  to  the  tick  marks  on  the  range  axis. 
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DEPTH  AVERAGED  LOSS 


Range  (km) 

SAFARI-FIP  case  3. 

Fig.  11.  Results  for  SAFARI-FIP  esse  3:  (s)  depth-averaged  transmission 
loss  as  function  of  range;  (b)  contour  plot  of  transmission  loss  as  function  of 
range  and  depth. 
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(a)  TRANSMISSION  LOSS 


(b)  TRANSMISSION  LOSS 


Fig.  12.  Results  for  SAFARI-PIP  case  3:  (a)  trans¬ 
mission  loss  vs  depth  at  2  km  range;  (b)  transmis¬ 
sion  loss  vs  depth  at  4  km  range. 
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■  7.3.4.  SAFARI-FIP  case  <?;  Long-range  propagation 

In  this  example  it  is  demonstrated  how  long-range  propagation  problems  are  treated 
by  properly  truncating  the  integration  interval.  The  problem  considered  is  exactly 
the  same  as  in  case  2,  but  now  the  transmission  loss  should  be  calculated  out  to  a 
range  of  50  km.  If  the  wavenumber  interval  was  left  unchanged,  several  thousand 
sampling  points  would  be  required.  It  is  well  known,  however,  that  except  for 
very  short  ranges  -  typically  less  than  a  few  water  depths  -  the  field  will  be  entirely 
dominated  by  the  two  propagating  normal  modes,  clearly  showing  up  as  sharp  peaks 
in  the  integrand  plot  of  Fig.  10a.  The  long-range  propagation  loss  can  therefore  be 
calculated  by  including  only  the  discrete  part  of  the  spectrum  containing  the  two 
modes,  significantly  reducing  the  required  sampling.  This  is  done  by  creating  the 
data  file  given  in  Table  14. 

The  resulting  plots  of  the  integration  kernel  and  the  transmission  loss  are  shown  in 
Fig.  13.  As  can  be  observed,  the  wavenumber  interval  has  been  properly  truncated. 
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Table  14 

SAFARI-FIP  cue  4  data  file 


Data  file 

Description 

Notes 

sm&l-FXP  case  4 

title 

»  I  T  J 

options 

90  0 

frequency  and  contour  offset 

6 

number  of  layers 

0  0  0  0  0  0  0 

layer  1  (upper  halfspace,  vacuum) 

0  1600  -1480  0010 

layer  2  (water  layer) 

30  1480  -1490  0  0  1  C 

layer  3  (water  layer) 

100  1600  400  0.2  0.6  1.8  0 

layer  4  (silt  sediment) 

120  1800  600  0.1  0.2  2.0  0 

layer  5  (sand  sub-bottom) 

60 

source  depth  50  m 

100  100  1  1 

receiver  depth  100  m 

1300  1760 

phase  velocities  CHIN  and  CMAX 

l 

612  46  466 

wavenumber  sampling  parameters 

2 

0.0  50  20  10 

range  axis  parameters 

40  120  12  20 

transmission-loss  axis  parameters 

1  A  proper  value  of  CNIX  is  found  by  observing  the  full  integrand  in  Fig.  10a.  An 
obvious  point  to  truncate  is  where  the  amplitude  vanishes  at  the  wavenumber  0.1 1  m"1, 
corresponding  to  a  phase  velocity  of  1712  m/s.  The  actual  value  of  CMS  has  been 
chosen  a  little  larger,  as  1750  m/s.  The  value  of  CHIN  is  less  critical,  but  1300  m/s  has 
been  chosen  to  yield  a  reasonable  range  sampling. 

1  Since  CHiX  was  not  chosen  equal  to  the  actual  value  corresponding  to  the  wavenumber 
where  the  integrand  amplitude  was  zero,  IC1  is  given  a  value  which  will  start  calcu¬ 
lation  at  this  particular  point.  Also,  since  the  integrand  amplitude  is  known  to  be 
insignificant  at  the  high  wavenumber  end  of  the  selected  interval,  IC2  is  set  to  a  value 
smaller  than  NV. 
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(a)  INTEGRAND 


(b)  TRANSMISSION  LOSS 


Fig.  13.  Results  for  SAFARI-FIP  case  4:  (a)  trun¬ 
cated  integration  kernel  as  function  of  horizontal 
wavenumber;  (b)  calculated  transmission  loss  as  func¬ 
tion  of  range. 
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■  7.5.5.  SAFARI-FIP  cue  5:  Beam  propagation 

In  this  last  example  it  is  demonstrated  how  SAFARI-FIP  is  applied  to  a  beam  prop¬ 
agation  problem.  A  linear  vertical  array,  placed  at  mid- water  depth  in  the  environ¬ 
ment  treated  in  case  1,  is  generating  a  1000  Hs  gaussian  beam  impinging  on  the 
bottom  at  a  nominal  gracing  angle  of  25°.  The  task  is  to  investigate  the  reflection 
and  transmission  characteristics  of  this  beam  by  generating  a  contour  plot  of  the 
sound-pressure  field  in  depth  and  range.  The  data  file  is  set  up  in  Table  15. 

The  contour  plot  produced  by  means  of  the  DISSPLA  plot  package  is  shown  in 
Fig.  14a,  whereas  Fig.  14b  shows  the  corresponding  UNIRAS  contour  plot  with  dark 
shading  indicating  high  levels.  Since  the  angle  of  incidence  of  the  beam  is  between  the 
critical  angles  of  the  silt  layer  and  the  sub-bottom,  the  incident  beam  is  both  reflected 
and  transmitted  at  the  water /sediment  interface  whereas  the  transmitted  beam  is 
totally  reflected  at  the  sediment/subbottom  interface.  The  resulting  complex  beam- 
splitting  reflectivity  pattern  is  easily  observed.  More  beam  examples  can  be  found 
in  [21,22,26,27]. 
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Table  15 

SAFARI-FIP  case  5  data  file 


Data  file 

Description 

Notes 

Smil-FXF  case  5 

title 

P  N  C  I. 

options 

1 

1000 

frequency  1  kHs 

4 

number  of  layers 

2 

65  1600  00010 

layer  1  (upper  halfspace,  water) 

2 

65  1500  00010 

layer  2  (isovelocity  water  layer) 

2 

100  1600  400  0.2  0.5  1.6  0 

layer  3  (silt  sediment) 

120  1800  600  0.1  0.2  2.0  0 

layer  4  (sand  sub-bottom) 

50  41  0.76  26.0  4  100 

source  array  parameters 

3 

60  126  61  62 

receiver  depths  50  to  125  m 

1600  6000 

phase  velocities  CK1N  and  CH1X 

4 

2046  1400  1960 

wavenumber  sampling  parameters 

5 

0.0  0.3  20  0.06 

range  axis  for  contour  plot 

50  126  12  26 

depth  axis  for  contour  plot 

24  64  6 

contour  levels  in  dB 

1  Option  P  indicates  that  the  sources  are  line  sources  rather  than  point  sources.  The 
calculations  are  therefore  performed  in  plane  geometry,  and  the  calculated  parameter 
is  the  normal  stress,  option  N.  Option  C  generates  a  contour  plot  of  the  field  in  depth 
and  range.  Finally,  option  L  indicates  that  the  field  is  generated  by  a  vertical  linear 
array. 

1  Since  only  the  first  bottom  bounce  of  the  beam  is  of  interest,  the  sea  surface  is  remove.! 
by  replacing  the  water  column  by  an  infinite  halfspace.  Note  that  a  dummy  interface 
has  been  introduced  just  below  the  lowermost  source  in  the  array.  This  is  a  ‘trick’, 
reducing  the  computation  time  in  cases  where  many  sources  and  receivers  are  present. 
This  is  due  to  the  fact  that  the  homogeneous  solution  within  a  layer  has  to  be  super¬ 
imposed  with  the  direct  source  field  -  an  operation  which  must  be  performed  for  every 
aource/receiver  combination  within  each  layer.  In  the  present  cas->  the  introduction  of 
the  dummy  interface  reduces  the  number  of  receivers  present  in  the  source  layer  to  the 
uppermost  10,  and  the  savings  in  computational  effort  more  than  compensates  for  the 
additional  computation  involved  in  adding  an  extra  interface.  It  is  left  as  an  exercise 
to  the  user  to  demonstrate  this. 

3  Since  option  L  was  chosen,  all  7  source  parameters  have  to  be  specified.  The  mid-point 
of  the  array  is  50  m,  i.e.  in  the  middle  of  the  wate-  column.  The  array  has  41  elements 
with  a  spacing  of  0.75  m,  i.e.  half  a  wavelength.  The  total  array  length  is  therefore 
30  m.  The  steerir  j  angle  is  25*  downward  with  respect  to  horisontal.  The  array  is  of 
type  4,  indicating  a  Gaussian,  non-focusing  shading.  The  focal  depth,  specified  to  the 
depth  of  the  sea-bed,  is  therefore  dummy  in  this  case 

*  Although  the  wavenumber  spectrum  of  the  parallel  beam  will  be  relatively  narrow, 
almost  the  whole  spectrum  of  propagation  angles  in  the  water  has  been  included  in 
order  to  get  a  reasonable  range  sampling. 

1  The  wavenumber  sampling  parameters  have  been  chosen  such  that  only  the  non¬ 
vanishing  part  of  the  spectrum  is  included.  This  can  be  shown  by  adding  option 
I  to  obtain  plots  of  the  integration  kernel. 
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8.  Running  SAFARI-FIPP 

The  SAFARI-FIPP  module  calculates  the  depth-dependent  Sreen’s  function  for  a 
selected  number  of  frequencies  and  determines  the  pulse  response  at  a  given  re¬ 
ceiver  position  by  evaluating  first  the  wavenumber  integral,  Sq.  (13),  and  then  the 
frequency  integral,  Eq.  (8).  As  was  the  case  for  SAFARI-FIP,  both  stresses  and  par¬ 
ticle  velocities  can  be  determined,  and  the  field  may  be  produced  by  either  point 
or  line  sources.  By  arranging  the  sources  in  a  vertical  phased  array,  pulsed  beam 
propagation  can  be  analysed. 

The  time  response  is  of  course  dependent  on  the  pulse  shape  at  the  source.  In 
SAFARI-FIPP  the  user  can  either  select  one  of  the  five  internal  pulse  shapes  or  create 
an  external  file  containing  any  desired  pulse  shape.  In  either  case  the  source  pulse 
is  defined  as  the  pressure  pulse  produced  at  a  distance  of  1  m  from  the  source  (for 
solids  the  negative  of  the  normal  stress  1  m  below  the  source).  The  time  dependence 
is  therefore  not  directly  that  of  the  forcing  term  in  the  potential  wave  equation, 
Eq.  (6).  This  is  important  to  note  when  comparing  the  SAFARI-FIPP  solutions  with 
those  obtained  by  other  codes. 

The  pulse  response  output  is  available  as  either  individual  pulse  plots  for  single 
receivers  or  as  stacked  plots,  where  the  stacking  can  be  performed  in  either  range  or 
depth.  For  the  individual  plots  the  time  series  are  produced  in  true  units,  i.e.  Pa  for 
stresses  end  m/s  for  particle  velocities,  again  assuming  that  the  source  pulse  shape 
is  given  in  Pa.  For  the  stacked  plots  the  individual  traces  are  scaled  according  to 
certain  rules  specified  in  the  following. 

All  inputs  to  SAFARI-FIPP  are  read  from  the  file  currently  assigned  to  the  logical 
file  FOROOl.  Before  running  the  program  the  user  has  to  assign  the  file  containing 
the  input  data  to  this  logical  name.  The  most  convenient  approach  is  to  incluue  this 
assignment  in  a  general  command  file  which  also  assigns  file  names  to  the  logical 
names  of  the  output  files,  as  described  in  Subsect.  5.6. 

As  was  the  case  for  the  single-frequency  module,  successful  use  of  SAFARI-FIPP  is 
highly  dependent  on  the  parameters  specified,  in  particular  those  related  to  the 
truncation  and  discretization  for  carrying  out  numerical  integrations.  It  is  obvious 
that  the  additional  frequency  integration  performed  in  SAFARI-FIPP  makes  this 
module  even  more  difficult  to  use  than  SAFARI-FIP.  It  is  therefore  important  that 
the  user  be  extremely  confident  in  running  SAFARI-FIP  before  trying  to  use  the  pulse 
model.  In  the  following  the  preparation  of  input  files  is  first  discussed  in  detail  and 
then  the  numerical  considerations  particular  to  the  pulse  model  are  addressed. 
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Table  16 

Parameters  of  SAFARI-FIPF  input  files:  Calculation  and  environmental  parameters 


Block 

Parameter 

Unite 

Limits 

I 

TITLE  : 

title  of  run 

- 

<  80  char. 

II 

optl  opt2 

output  options 

- 

<  40  char. 

III 

FREQ: 

source  centre  frequency 

Hs 

>  0 

COFF: 

integration  contour  offset 

dB/A 

COFF  >  0 

IV 

ML: 

number  of  layers,  incl.  halfspaces 

- 

ML  >  2 

D: 

depth  of  interface 

m 

- 

CC: 

compressional  speed 

m/s 

CC  >  0 

CS: 

shear  speed 

m/s 

- 

AC: 

compressional  attenuation 

dB/A 

AC  >  0 

AS: 

shear  attenuation 

dB/A 

AS  >0 

RO: 

density 

g/cmJ 

RO  >  0 

RG: 

rms  value  of  interface  roughness 

m 

- 

CL: 

correlation  length  of  roughness 

m 

CL  >  0 

V 

SD: 

source  depth  (mean  for  array) 

m 

— 

NS: 

number  of  sources  in  array 

- 

MS  >  0 

DS: 

vertical  source  spacing 

m 

DS  >  0 

AN: 

gracing  angle  of  beam 

deg 

- 

I*: 

array  type 

- 

1  <  IA  <  5 

FD: 

focal  depth  of  beam 

m 

FD  ^  SD 

VI 

RD1: 

depth  of  first  receiver 

m 

— 

RD2: 

depth  of  last  receiver 

m 

RD2  >  RD1 

MR: 

number  of  receivers 

- 

NR  >  0 

VII 

CHIN: 

minimum  phase  velocity 

m/s 

CHIN  >  0 

CM  AX: 

maximum  phase  velocity 

m/s 

- 

VIII 

MV: 

number  of  wavenumber  samples 

— 

MV  >  1 

IC1: 

first  sampling  point 

- 

IC1  >  1 

IC2: 

last  sampling  point 

- 

IC2  <  MV 

IINC: 

integrand  plot  increment 

IINC  >  0 

IX 

NT: 

number  of  time  samples 

MT=  2U 

FI: 

low  frequency  limit 

Hs 

FI  >  0 

F2: 

high  frequency  limit 

Hs 

FS  >F1 

DT: 

time  sampling  increment 

s 

DT  >  0 

RO: 

first  receiver  range 

km 

- 

DR: 

receiver  range  increment 

km 

DR  ^  0 

MEAN: 

number  of  ranges 

- 

MEAN  >  0 
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8.1.  INPUT  FILES  FOR  SAFARI-FIPP 

The  input  data  are  structured  in  13  blocks.  The  first  nine,  described  in  Table  16, 
specify  the  title,  options,  frequency,  environmental  parameters,  source  and  receiver 
geometry  and  finally  the  wavenumber  and  frequency  integration  parameters.  The 
last  four  blocks,  described  in  Table  17.  contain  axis  specifications  for  the  graphical 
output.  Some  of  these  blocks  should  always  be  included  and  others  only  if  certain 
options  have  been  specified.  The  single  blocks  and  parameters  are  described  in  detail 
in  the  following. 


Table  17 

Parameters  of  SAFARI-FIPP  input  files:  Plot  parameters 


Block 

Parameter 

Units 

Limits 

X1 

CRED: 

reduction  velocity 

m/s 

CRED  >  0 

THIN: 

start  of  time  window 

8 

- 

THAI: 

end  of  time  window 

8 

- 

TLEN: 

length  of  time  aixs 

cm 

TLEN  >  0 

TINC: 

distance  between  tick  marks 

8 

TINC  >  0 

XIJ 

SPLO: 

lower  limit  of  stacked  plots 

m  or  km 

— 

SPUP: 

upper  limit  of  stacked  plots 

m  or  km 

- 

SPLN: 

length  of  stacking  axis 

cm 

SPLN  >  0 

SPIN: 

distance  between  tick  marks 

m  or  km 

SPIN  >  0 

XII3 

NMOD: 

number  of  modes 

- 

NMOD  >  1 

FNIN: 

min.  freq.  dispersion  curves 

He 

FMIN  >  0 

FMAX: 

max.  freq.  dispersion  curves 

Hi 

FMAX  >  0 

FLEN: 

lenght  of  frequency  axis 

cm 

FLEN  >  0 

FINC: 

frequency  tick  mark  spacing 

He 

FINC  >  0 

GVLO: 

min.  phase/group  velocity 

ra/s 

GVLO  >  0 

GVUP: 

max.  phase/group  velocity 

m/s 

GVOP  >  0 

GVLN: 

lenght  of  velocity  axis 

cm 

GVLN  >  0 

GVIN: 

velocity  tick  mark  spacing 

m/s 

GVIN  >  0 

XIII4 

VLEF: 

wave  speed  at  left  border 

m/s 

— 

VRIG: 

wave  speed  at  right  border 

m/s 

- 

VLEN: 

length  of  wave  speed  axis 

cm 

VLEN  >  0 

UNC: 

wave  speed  tick  mark  distance 

m/s 

VINC  >  0 

DV’  : 

depth  at  upper  border 

m 

- 

DVLJ: 

depth  at  lower  border 

m 

- 

DVLN: 

length  of  depth  axis 

cm 

DVLN  >  0 

DVIN: 

depth  axis  tick  mark  interval 

m 

DVIN  >  0 

1  Only  for  option  &  or  NRAN  >  0.  3  Only  for  options  S  and  D.  3  Only  for  option  G. 

4  Only  for  option  Z. 
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I 

i 


BLOCK  XT  -  OPTIONS 

Lines:  i 

Requited:  Always 

Syntax:  opt 1  opt 2  opt3. . . 

JFbnnct;  Free.  Maximum  40  alphanumeric  characters. 

Description:  Series  of  alphanumeric  characters  specifying  the  desired 

options.  The  following  are  implemented  at  present: 


B:  If  this  option  is  specified  the  slowness  interval 
(or  phase  velocity  interval)  will  be  held  constant 
at  all  frequencies.  As  a  default  procedure, 
the  wavenumber  interval  will  be  determined 
from  the  specified  phase  velocities  and  the 
maximum  frequency  and  will  then  be  held 
constant  at  all  lower  frequencies.  The  B  option 
should  be  used  when  only  a  certain  part  of  the 
Singular  spectrum  is  considered,  e.g.  the  discrete 
spectrum  containing  normal  modes.  On  the  other 
hand,  When  evanescent  waves  are  considered, 
e.g.  seismic  interface  waves,  the  wavenumber 
interval  should  be  held  constant  in  order  to 
Include  the  whole  exponential  ‘tail’  at  the  low 
•  frequencies;  thus  in  these  cases  the  B  option 
should  he  omitted.  It  is  always  a  good  habit 
to  check  the  integrand  at  a  selected  number  of 
frequencies  (UHC  >  0) 

D:  Fbr  each  selected  range  the  pulses  will  be 
plotted  in  a  depth-stacked  format.  This  is  the 
VSP  (vertical  seismic  profiling)  type  of  output 
used  in  seismology.  No  geometrical  spreading 
compensation  is  performed  in  the  depth-stacked 
format,  but  a  common  sealing  factor  determined 
for  the  upper  trace  is  applied.  As  for  option  S, 
the  scaling  is  performed  individually  for  each 
stacked  plot;  the  amplitudes  for  different  ranges 
therefore  cannot  be  directly  compared.  The 
depth-stacked  plots,  however,  can  be  rescaled  by 
running  FIPPLOT  with  option  RES. 

F  t  The  Filon  integration  scheme  is  used  for  the 
wavenumber  integration  instead  of  the  standard 
scheme  based  on  the  trapezoidal  rule. 
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BIOOC  XT  -  (Continued) 


Qm&ptioa: 


V-  \v  \ :  v : 


Diap*  .on  curve*  are  plotted  for  a  selected 
number  of  interface  modes  and  normal  modes. 
The  modes  sure  not  found  directly  by  searching  for 
singularities  but  indirectly  by  localising  maxima 
in  the  depth- dependent  Green’s  function.  The 
user  should  therefore  be  sure  that  the  mode  of 
interest  is  present  at  the  selected  depth  at  all 
frequencies.  Otherwise  a  ‘false’  mode  may  be 
found,  leading  to  discontmous  dispersion  curves. 
Horisonta!  particle  velocity  is  calculates 
Complex  wavenumber  integration  :ontour  is 
applied.  The  contour  is  shifted  into  the  upper 
halfplane  by  an  offset  controlled  by  the  input 
parameter  COFF  (Block  ID). 

Sources  are  arranged  in  a  vertical  line  array. 
Normal  stress  <r„  (=  -p  in  fluids)  is  calculated. 
Plane  geometry.  The  sources  will  be  line  source* 
instead  of  point  sources  as  used  in  the  default 
cylindrical  geometry.: 

Plpt  of  source  pulse  at  1  m  distance  from  a  single 
source  will  be  produced  together  with  a  plot  of 
the  associated  frequency  spectrum.  The  chosen 
source  pulse  is  filtered  to  exclude  frequencies 
outside  the  interval  considered  in  the  calculation 
(Block  IX'  ). 

For  each  selected  receiver  depth  the  pulses  are 
plotted  in  a  range-  stacked  format.  The  vortical 
axis  of  the  plots  represent*  the  range,  and  each 
receiver  tr  ee  is  placed  with  the  sero-liae  at  the 
corresponding  range  value.  In  order  to  remove 
the  effect  of  spherical  spreading,  the  amplitudes 
are  first  multiplied  by  the  range  of  the  actual 
receiver,  and  then  a  common  scaling  factor  is 
applied  to  all  traces.  This  factor  is  determined 
such  that  the  trace  for  the  first  receiver  range 
does  not  overlap  the  next  trace.  This  scaling  is 
performed  individually  for  each  receiver  depth; 
thus  the  amplitudes  for  different  receivers  cannot 
be  directly  compared.  The  scaling  may  be 
changed  by  running  the  plot  program  FIPPLOT 
with  the  R£S  option  as  described  in  Appendix  B. 
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SaOLANTCEN  SR- 113 


BLOCK  U  -  (Continued) 

Detcriptioa: 

?:  Vertical  particle  velocity  ia  calculated. 

X:  If  sources  are  present  in  a  solid  medium,  these 
will  be  considered  vertical  point  forces  (line  forces 
in  plane  geometry)  rather  than  omnidirectional 
sources. 

Z:  Plot  of  velocity  profile  is  generated, 
a:  Digit  determining  the  desired  source  pulse, 

defined  as  the  pressure  or  normal  stress  produced 
at  1  m  distance  from  the  source.  There  are 
five  internal  pulse  shapes,  shown  together  with 
their  normalised  frequency  spectra  in  Figs.  15 
to  19.  The  user  can,  however,  create  any  desired 
pulse  shape  in  an  external  file  and  read  it  in  by 
specifying  n  =  0.  The  internal  pulse  shapes  in 
pressure  (negative  of  normal  stress)  are  functions 
of  time  t  and  the  centre  angular  frequency  b-c  ~ 

2 *fQ  as  follows: 

(1]  This  is  the  default  source  pulse,  defined  by  the 
function 

/(t)  «  0.75  ~  cosu>ct  +  0.25  cos  2u>ct, 

0<i<r«l//e. 

[2]  This  source  pulse  has  the  shape 

s  ; 

/(<)  *  ^AiPton^irit/T), 

0  <  t  <  T  =  1.55/A. ' 

where  A\  =  0.48829,  A*  -  -0.14128  and 
—  0.01158,  ensuring  vanishing  first  and  second 
derivatives  at  t  =  0  and  t  =  T. 
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SLOCK  XT  - 

(Continued) 

Description- 

[3]  A  single  sine  period 

f(t)  =  sinu>et,  0  <  t  <  T  =  l//c. 

[4]  This  option  creates  a  source  pulse  with  a  time 
dependence  consisting  of  4  sine  periods,  weighted 
with  a  Hanning  window, 

f(t)  =  jr  sino)et  (l  -  cos  \u>ct)  . 

0  <  t  <  r  =  4/  fc. 

[5]  The  last  internal  source  pulse  is  the  one  often 
used  in  seismology  [14], 

f{t)  s=  sinu>e<  ~  !-sin2u\.t, 

0  <  t  <  T  «  l//c. 

If  option  [0)  is  specified  the  source  pulse  is  read 
the  file  currently  assigned  to  the  logical  unit 
j  66.  This  fue  should  be  formatted  and  contain 
;  the  amplitude  sampled  with  the  same  interval  as 
specified  by  the  parameter  DT  (Block  IX),  and 
with  one  sample  per  record.  The  duration  of 
the  source  pulse  is  determined  automatically  by 
detection  of  the  end-oMile.  : 

Pressure  (  Pa  ) 
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Saclantcen  sr- ns 


(a)  SQ.URCE  -EULSE 


(b)  SOURCE.  SEECIRUM 


Fig.  IS.  SAFARI-FIPP  source  poise  no-  1  and  its 
frequency  spectrum. 
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(a)  SOURCE  PULSE 


(b)  SOURCE  SPECTRUM 


Fig.  Jo.  SAFARI-  CIPP  source  pulse  no.  «  and  its 
frequencv  spectrui.i. 
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SACLANTC2N  SR.  US 


(a)  SOURCE  PULSE 


(b)  SOURCE-SPECTRUM 


0.0  1.0  <0  3.0  4.0  5.0 

Dimensionless  frequency  f/f« 

Fig.  17.  SAFARI-FIPP  source  pulse  no.  3  end  its 
frequency  spectrum. 
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(a)  SOURCE..  P-ULSE 


(b)  SOURCE  SPECTRUM 


18.  SAFARI-FIPP  aource  pulse  no.  4  and  its 
frequency  specttam. 


8  -  Running  SAFARl-FIPP 


(a)  SPURGE  PULSE 


(b)  SOURCE  SPECTRUM 


Fig.  Ji.  SAFARl-FIPP  source  pulse  no. 
frequency  spectrum. 
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BLOCK  m  * 

FREQUENCY 

line*: 

1 

Required: 

Always 

Syntax; 

FREQ,  COFF 

format: 

Free 

Detcriptioo: 

FREQ:  Centre  frequency  fe  of  source  spectrum  in  Hs. 

COFF:  Wavenumber  integration  contour  offset.  Tb  be 
specified  in  dB/A,  where  A  is  the  wavelength  at 
the  source  depth  SD.  Since  only  the  horizontal 
part  of  the  integration  contour  is  considered,  this 
parameter  should  not  be  chosen  so  large  that  the 
amplitudes  at  the  ends  of  the  integration  interval 

s  . 

become  significant.  In  lossless  cases  too  small 

values  will  give  sampling  problems  at  the  normal 
modes  and  other  singularities.  For  intermediate 
values,  the  result  is  independent  of  the  choice 
of  COFF,  but  it  is  recommended  to  use  a  value 

which  at  the  centre  frequency  fe  is  the  same  as 
that  used  by  SAFARI*F!P  for  the  same  number  of 
sampling  points,  i.e. 

COPP  =  60cc(*,)(e-}n  ~  c-U)/NW. 

'S^V'-'VV’''  ■'  - 

!•:  V: :  .  .  •••  : 

!,W'  .  $  ’  ■<  ‘  •  •  ;  ?  s  .  ; 

This  value  is  the  default  which  is  applied  if  option 
J  is  specified  and  COFF  =  0. 
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SACLANTCBN  SR-113 


BLOCK  XV  -  ENVIRONMENTAL  DATA 

Lines:  >3 

Required;  Always 

Syntax, 

NL 

D(l)  CC(1)  CS(1)  AC<1)  AS(l)  R0(1)  RG(1)  CL(1) 
0(2)  CC(2)  CS(2)  AC (2)  AS(2)  R0(2)  RG(2)  CL(2) 


Description: 


tit: 


DO: 


CCO: 


CS(): 


ACO: 


AS(> 


Number  of  layers,  including  the  upper  and  lower 
half0paces.  These  should  always  be  included,  even 
in  cases  where  they  are  vacuum. 

Depth  z  in  ro  of  upper  boundary  of  layer  or 
halfspace.  The  reference  depth  can  be  choosen 
arbitrarily,  and  D()  is  allowed  to  be  negative. 

For  layer  no.  1,  i.e.  the  upper  halfspace,  this 
parameter  Is  dummy. 

Velocity  c«  of  compressional  waves  in  m/s. 

If  specified  as  0.0,  the  layer  or  halfspace  is  a 
vacuum. 

Velocity  c,  of  shear  waves  in  m/s.  In  order  to  he 
physically  meaningfal,  it  is  required  that  t,  < 
v/O  JSCt.  If  specified  as  0.0,  the  layer  or  halfspace 
is  fluid.  K<0  <  0,  it  represents  the  compressional 
velocity  at  the  bottom  of  the  actual  layer,  which 
is  treated  as  fluid  with  l/c(*)*  varying  linearly 
with  depth. 

Attenuation  >  of  compressional  waves  in  dB/A. 

If  the  layer  is  fluid,  and  AC<)  is  specified  as  0.0, 
then  an  empirical  water  attenuation  is  used. 
Attenuation  y,  of  shear  waves  in  dB/A.  For  the 
attenuations  to  be  physically  meaningful,  it  is 
required  that 

It 

7c 
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BLOCK  tV  -  (Continu'd) 

Dmcsiptioa: 

R0() :  Density  p  ing/cm9. 

ROC)  i  RMS  roughness  of  interface  in  m.  RG(1)  is 

dummy.  If  RG  >  0  the  Kirchhoff  approximation 
is  used  forthe  scattering  loss,  i.e»  CL  -4  00. 

This  is  computationally  faster  than  doing  the 
full  scattering  theory,  but  it  is  oMy  applicable  for 
small  roughnesses  and  large  correlation  lengths, 
Tb  invoke  the  full  non-Kirchhoff  scattering  RG 
must  be  specified  equal  to  the  negative  rms 
roughness,  i.e.  RG  <  0. 

CL() :  Correlation  length  of  roughness  in  m.  This 
parameter  is  only  required  if  non- Kirchhoff 
scattering  has  been  selected,  i.e.  RG  <  0. 
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BLOCK :  V  - 

SOURCE  DATA 

Lines: 

1 

Required: 

Always 

Syntax: 

Option  L:  SD  NS  DS  AN  IA  FD 
else:  SD 

Format: 

Free 

Description: 

SD:  Source  depth  z*  in  m.  If  option  L  has  been 

specified,  SD  defines  the  mid  point  of  the  vertical 
source  array. 

NS:  Number  of  sources  in  the  array. 

DS :  Source  spacing  in  m. 

AN:  Specifies  the  nominal  grafting  angle  of  the 

generated  beam  in  degrees.  AN  >  0  corresponds 

fit hniUNNt 

XA:  Array  type 

[1]  Rectangular  weighted  array 
[2}  Hanningweightedarray 

[3]  Hanning  weighted  focused  array 

[4]  Gaussian  weighted  array 

[5]  Gaussian  weighted  focused  array 

FD :  Focal  depth  in  ra  for  an  array  of  type  3  or  5. 
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BIOCK  VI  - 

RECmmPATA 

lines; 

1 

Requited: 

Always 

Syntax: 

RD1  RD2  NR 

Fbrraat. 

>:  Free 

Description: 

WU:  Receiver  depth  in  m.  If  the  calculations  are  to  be 

performed  for  more  than  one  receiver,  i.o.  NR  >  1  , 
then  RD1  specifies  the  uppermost’ receiver. 

RD2:  If  NR  >  1  this  parameter  determines  the  aepth 
of  the  lowermost  receiver.  If  NR  «  1,  then  the 

parameter  twz  Qummy. 

NR:  Number  of  receiver  depths  considered. 
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BWQK  vn  -  PBASB  VELOCITIES 

1 

Always 
CHIN  CNAI 
free 

CHIU:  Minimum  phene  velocity  cmin  in  m/s.  Determines 
the  upper  limit  of  the  truncated  horizontal 
wavenumber  space: 

i’nni  "  2x/ / Cmin 

where  /  is  the  actual  frequency  if  option  B  was 
selected;  otherwise  /  is  the  maximum  frequency 
F2.  In  this  last  case  will  therefore  be 
'  independent  of  frequency, 

CMAlt  Maximum  phase  velocity  cmM  in  m/s. 

Determines  the  lower  limit  of  the  truncated 
horizontal  wavenumber  space: 

hmlft  ~  2tr//fnuM: 


where  /  is  the  actual  frequency  if  option  B  was 
selected;  otherwise  /  is  the  maxinr’m  frequency 
F2.  In  this  last  case  will  therefore  be 
independent  of  frequency.  In  plane  geometry 
(option  P)  QUX  may  be  specified  as  negative,  and 
the  negative  wavenumber  spectrum  will  then  be 
included  with  yielding  a  correct 

solution  also  at  aero  range.  Contour  offset  (option 
Jf)  is  not  allowed  in  this  case. 


Idtmi 
Aoqttfred; 
Syntax;  ■ 
Fbdnat: 
Description; 
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MQCZ  m 


wavenumber  sampling 


Number  of  sampling  point*  in  wavenumber  space. 
Since  the  Fait  Field  technique  is  not  used  for  the 
wavenumber  integration,  JIW  may  be  any  positive 
integer.  By  specifying  NW  a*  1  the  propagation 
of  a  single  plane  or  conical  wave  component 
may  be  analysed.  Apart  from  this  trivial  case, 
the  sampling  points  are  placed  equidistantly  in 
the  truncated  wavenumber  space  determined 
by  CHIN  and  CM  AX  ,  If  the  default  trapesoidal 
rule  integration  is  applied,  it  is  recommended 
to  choose  MW  so  large  that  the  wavenumber 
sampling  Akc  at  the  centre  frequency  satisfies  the 
inequality 

<  i* 


where  rmtt  is  the  maximum  range  considered  in 
the  calculations.  For  the  Filon  scheme  a  smaller 
number  of  sampling  points  may  be  chosen,  hut 
this  depends  on  how  smooth  the  integration 
kernel  is.  In  any  case  the  sampling  should  be 
checked  by  convergence  tests. 

Number  of  the  first  sampling  point,  at  which 
the  calculation  is  to  he  performed.  IfXCl  >  1, 
then  the  Henkel  transform!*  saroed  for  aampfia* 

points  1,2  ...  IC1-1,  and  the  discontinuity  is 
smoothed. 

Number  of  the  last  sampling  point  where  the 
calculation  is  tp  be  performed.  If  IC2  <  Htf  , 
then  the  Hankel  transform  is  seroed  for  sampling 
points  IC2  4*1, ...  MW,  and  the  diseontinuitv  it 
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BtQCtC  m  - 

TIME,  FREQUENCY  AND  RANGE  PARAMETERS 

Linear 

1 

Requited: 

Always 

Syntax: 

Hi  PI 

F2  DT  RO  DR  NRAN 

Fbrmat: 

free 

Description; 

MT: 

Number  of  sampling  points  in  time  domain. 

Should  be  an  integer  power  of  2.  To  avoid 
aliasing  it  is  extremely  important,  that  the  time 

Pis 

window  (NT*DT)  is  chosen  large  enough  to  contain 
the  entire  time  signal  at  all  the  selected  ranges. 
Lower  limit  of  frequency  band  considered.  Should 

F* :  • 

be  given  in  Hs. 

Upper  limit  of  frequency  band  considered.  Should 

be  given  in  Hs.  If  F2  is  specified  as  having  a 

value  greater  than  the  Nyquist  frequency  0.5/DT, 

it  will  be  replaced  by  the  Nyquist  frequency 

in  the  calculations.  It  should  be  noted  that  a 
simple  rectangular  weighting  is  applied  to  the 
frequency  window  defined  by  PI  and  F2  .  To 
avoid  'ringing’  in  the  pulse  responses  due  to  the 
abrupt  truncation,  it  is  therefore  required  that 

PI  and  F2  be  chosen  where  the  source  frequency 
spectrum  has  a  minimum,  see  Figs.  15  to  19. 

FT: 

Interval  between  samples  in  time  domain.  Should 

;  '  ! 

be  given  in  seconds.  See  notes  above. 

-  ■'  ■ 

AO: 

Bnttge  in  km  for  first  receiver. 

FA: 

Spacing  in  km  between  receivers  (horisontal 

:f|| 

direction). 

XA1X: 

Number  of  ranges.  Pulse  plots  will  be  produced 

for  each  of  the  selected  field  parameters  at  each 

\  ,.  •  -N  S  N  ■*;  '  ;  *  .  • 

:  . 

receiver  range. 
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X  -  Tim  WINDOW  PARAMETERS 


Jfoguitod: 

Description: 


For  option  R  and  if  HRAR  >0 
CRM 

THJR  THAI  TL£N  TING 
Flree 

5hii  block  determines  the  time  window  applied  in  all 
calculations  and  plots  of  pulse  responses.  Bar  integrand 
plots  all  sealing  and  labeling  is  automatic. 

CRBDr  Reduction  velocity  cr  in  m/s,  The  pulse 

calculations  and  the  pulse  plotting  are  performed 
in  the  reduced  time  t«.  ss  t  —  r/c*..  Fo  avoid 
aliasing  the  reference  velocity  should  be  chosen 
<  larger  than  the  largest  group  velocity  in  the 
problem.  If  CRBD  is  specified  as  0.0,  the  reduced 
time  will  be  set  equdto  the  real  time,  t*  s»  tr 
THIS  s  Minimum  time  in  seconds.  Determinesthe 

starting  point  for  the  reduced  time  axis  in  pulse ; 
plots.  : 

THAI  :  Maximum  time  in  seconds.  Determines  the  end 
point  of  the  reduced  time  axis  in  pulse  plots. 

TLEN :  Length  of  time  axis  in  cm. 

TXSC ;  Time  interval  between  tick  marks  on  the  time 
axis  in  seconds. 
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BLOCK  XX  ~ 

STACKED  PLOT  PARAMETERS 

LUmr 

\ 

t 

Pot  options  D  and  S 

Syntax; 

SPLO  SPUP  SPLH  SPIN 

, Fbtraat: 

PVee 

A/OVviipilOO. 

SPLO:  Lower  limit  of  y-axis  in  stacked  pulse  plots.  If 

the  range-stacked  format  has  been  chosen  (option 
3)  then  SPLO  should  he  in  km.  For  the  depth- 
stacked  format  (option  0)  it  should  be  given  in  m. 

V.:  \  v:  •  ••  ...X 

SPLO  should  be  chosen  such  that  enough  space  is 
available  for  keeping  the  negative  amplitudes  of 

■  X  '  :'X  ■  '■  ;■  •  ... :  ••  •: 

the  lowermost  trace  within  the  plot  frame. 

SPOP:  Upper  limit  of  y-axia  in  stacked  pulse  plots.  Same 

- 

uiuv*  w  *»n*u  space  luuum  uv  avauaDiv  ior 

(die  trace  closest  to  SPUP 

SPUU  Length  in  cm  of  y-axis  of  stacked  pulse  plots. 

SPJI:  Interval  between  tick  marks  on  y-axis  of  stacked 

pulse  plots.  Same  units  as  SPLO  and  SPOP 
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BLOCK  XU  -  PLOT  DATA  /  DISPERSION  CURVES 


linw 

Requited: 

Syntax: 


Format.’  :.. 
Description: 


If  option  G  ha*  been  specified. 
tmoD 

FNIN  FKAX  FLEN  FINC 
GVLO  GVUP  GVLN  GVIN 
Free 


MMOD: 


FHIMt 

PHU; 

FLEN: 

FINC: 

GVlOi 


GVW: 

GVLN.' 

GVIN; 


Number  of  inode*  for  which  phase  and  group 
velocity  curves  are  to  be  plotted.  In  this  context 
the  term  ‘modes’  covers  not  only  the  normal 
modes,  but  all  significant  peaks  in  the  Hankel 
transform  integrand,  thus  also  including  seismic 
interface  modes.  The  modes  are  counted  in  the 
usual  way,  i.e.  from  low  phase  velocities  towards 
high  phase  velocities. 

Minimum  frequency  in  Hs.  Determines  the 
starting  point  of  the  x-axis, 

Maximum  frequency  in  Hs.  Determines  the  end 
point  of  the  e-axis. 

Length  of  frequency  axis  in  cm. 

Tick  mark  interval  for  frequency  axis. 

Minimum  velocity  in  m/s.  Determines  the 
starting  point  of  the  y-axis  for  dispersion  curve 
plots. 

Maximum  velocity  in  ra/s.  Determines  the  end 
point  of  the  y-axis. 

Length  in  cm  of  velocity  axis  for  dispersion  plots. 
Velocity  tick  mark  interval  in  m/s. 
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Saclantcen  SR- 113 


MQCK  xm  - 

PLOT  DATA  /  VELOCITY  PROFILE 

tines. 

2 

RnquimJ: 

If  option  2  wu  specified  in  Block  II. 

Syntnx: 

VUSF  TRIG  VLEM  VINC 

.format: 

DVUP  DTLO  DTLN  DTIN 

Rree 

Description: 

TLEP :  Left  border  value  of  velocity  in  m/s. 

TRIG  t  Right  border  value  of  velocity  in  m/s. 

TIER:  Length  of  velocity  axis  in  cm. 

VINC:  Distance  between  tick  marks  in  m/s. 

DVUP:  Depth  at  upper  border  of  plot  in  m. 

DTLO:  Depth  at  lower  border  of  plot  in  m. 

DTLN  :  Length  of  depth  axis  in  cm. 

DT*I:  Distance  in  m  between  tick  marks  on  lepth  axis. 
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8.2.  NUMERICAL  CONSIDERATIONS 

It  mutt  be  assumed  that  the  user  is  confident  with  the  numerical  problems  involved 
in  running  the  reflection  coefficient  and  transmission  loss  models  SAFARI- FIPR  and 
SAFARI-FIP  before  trying  to  use  the  much  more  complicated  pulse  code.  We  will 
therefore  at  this  point  concentrate  on  the  numerical  considerations  particular  to  the 
pulse  code.  For  guidelines  concerning  the  selection  of  environmental  models,  the 
user  is  referred  to  the  previous  sections. 

The  first  thing  the  user  has  to  select  is  a  proper  source  pulse,  which  can  either  be  one 
of  the  internal  pulses  or  one  created  by  the  user  and  placed  in  a  file  that  can  be  read 
by  SAFARI-FIPP.  After  selecting  the  centre  frequency  /e,  the  frequency  bandwidth, 
defined  by  the  parameters  FI  and  F2,  is  chosen  such  that  the  amplitude  of  the  source 
spectrum  is  small  at  the  truncation  points.  This  is  easily  done  by  inspection  of  the 
frequency  spectrum  of  the  source  pulse. 

Next,  the  time  sampling  DT  is  selected.  If  chosen  too  large,  the  pulse  plots  will 
become  very  ‘peaky’,  but  in  most  cases  10  sampling  points  per  period  at  the  centre 
frequency  fe  will  be  sufficient,  i.e.  DT  =  0.1  /  fe.  The  total  time  window  T  is  of  course 
dependent  on  the  total  number  of  time  samples  NT,  but  it  is  usually  convenient  to 
postpone  the  determination  of  this  parameter  until  later. 

Instead  we  select  at  this  point  the  parameters  controlling  the  wavenumber  sam¬ 
pling.  The  numerical  considerations  in  this  regard  are  identical  to  those  described 
for  SAFARI-FIP,  but  the  wide  frequency  band  is  a  complicating  factor.  First  an 
appropriate  slowness  interval,  defined  by  CHIN  and  CHAX,  should  be  selected,  which 
is  most  conveniently  done  in  the  following  way: 

1.  Select  three  characteristic  frequencies,  typically  the  centre  frequency  and  one 
at  each  end  of  the  frequency  interval. 

2.  Determine  the  necessary  wavenumber  interval  for  each  frequency  as  described 
for  SAFARI-FIP. 

3.  Now  determine  whether  the  wavenumber  interval  should  be  held  constant 
at  all  frequencies  (default)  or  should  vary  such  that  the  slowness  interval  is 
constant  (option  B)  and  select  values  of  CNIN  and  CNAX  accordingly. 

If  only  a  certain  spectral  part  (limited  and  frequency-independent  beamwidth)  is  be¬ 
ing  considered,  then  option  B  must  be  selected,  but  in  this  case  the  procedure  above 
is  used  to  determine  proper  truncation  points,  and,  aj  was  the  case  for  SAFARI-FIP, 
the  truncation  may  have  to  be  tapered  by  selecting  proper  values  of  ICi  and  IC2. 

The  number  of  wavenumber  sampling  points  is  a  parameter  which  is  even  more 
critical  in  the  present  case  than  for  transmission  loss  calculations  where  an  error  of 
a  fraction  of  a  dB  is  usually  acceptable.  This  is  due  to  the  fact  that  not  only  the 
amplitude  but  also  the  phase  of  each  frequency  component  has  to  be  accurate  in  order 
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to  yield  a  correct  pulse  response.  Although  a  convergence  test  could  be  performed 
by  using  SAFARI-FIP  at  a  few  selected  frequencies,  it  will  usually  be  necessary  to 
perform  the  convergence  test  on  the  full  pulse  problem.  For  the  inexperienced  user, 
however,  it  is  advisable  to  perform  a  convergence  test  for  SAFARI-FIP  at  least  at  the 
centre  frequency  to  obtain  reasonable  initial  estimates  of  the  number  of  wavenumber 
samples.  It  should  be  stressed  that  although  the  Filon  integration  scheme  (option  F) 
is  not  default  it  is  highly  recommended  due  to  the  increased  accuracy  and  reduced 
sampling  requirements. 

An  initial  convergence  test  for  the  wavenumber  sampling  is  most  conveniently  done 
together  with  the  one  necessary  for  determining  the  length  of  the  time  window, 
T  =  NT  *  DT.  This  is  done  in  the  following  way: 

1.  Select  a  narrow  frequency  window  arout.  1  the  centre  frequency,  typically  FI  = 
0.8 fe  and  F2  =  1.2/c. 

2.  Choose  CRED  and  NT  such  that  the  time  window  of  length  T  is  believed  to  con¬ 
tain  all  arrivals.  As  the  frequency  sampling  interval  is  inversely  proportional 
to  the  length  of  the  time  window,  A/  =  1/T,  NT  should  be  chosen  as  small 
as  possible  in  order  to  limit  the  calculation  time.  CRED  is  most  conveniently 
taken  to  be  the  highest  wave  speed  in  the  problem. 

3-  Select  pulse  type  4  which  has  the  narrowest  spectrum. 

4-  Run  SAFARI-FIPP  to  obtain  pulse  responses  at  all  ranges. 

5-  If  the  pulse  response  at  longer  ranges  become  ‘noisy’,  it  indicates  that  the 
wavenumber  sampling  is  insufficient .  Increase  the  sampling  and  re-run  until 
the  ‘noise’  disappears.  Remember  that  the  very  narrow  frequency  band  will 
yield  considerable  ‘ringing’  of  the  pulse,  which  should  not  be  confused  with 
the  undersampling  noise.  The  amount  of  ‘ringing’  of  the  source  pulse  can 
be  checked  by  selecting  option  R,  which  produces  a  plot  of  the  source  pulse 
filtered  by  the  rectangular  window  determined  by  FI  and  F2. 

6-  When  the  result  is  stable,  increase  the  time  window  by  changing  either  NT  or 
DT.  The  true  arrivals  will  be  invariant  to  the  change  in  time  window,  whereas 
the  arrivals  which  are  wrapped  around  will  change  position.  Repeat  until  the 
wrap-around  disappears. 

7-  The  pulse  type  is  now  changed  to  the  desired  one  and  the  frequency  window  is 
extended  accordingly,  usually  containing  the  entire  main  lobe  of  the  spectrum. 
The  full  response  is  then  calculated.  If  high-frequency  noise  again  appears  at 
the  long  ranges,  the  wavenumber  sampling  must  be  increased.  On  the  other 
hand  it  will  only  rarely  be  necessary  to  change  the  time  window  again. 

8-  As  usual:  When  the  result  is  stable,  it  is  the  final  result. 

It  is  clear  that  this  convergence  test  procedure  can  become  extremely  time  consum¬ 
ing  since  the  computation  time  is  proportional  to  both  the  number  of  frequencies 


-  128  - 


SaCLANTCBN  SR-11J 


8  -  Running  SAFARI-FIPP 


and  the  number  of  wavenumber  samples.  It  is  therefore  important  that  the  user  se¬ 
lects  reasonably  low  initial  values,  which  again  requires  that  he  is  already  confident 
with  the  wavenumber  sampling  concepts  of  the  single  frequency  SAFARI-FIP  mod¬ 
ule.  Further,  it  is  extremely  important  that  the  user  has  a  substantial  knowledge 
of  time/frequency  analysis  by  means  of  FFTs  in  order  to  be  able  to  properly  select 
the  time/frequency  parameters  and  to  pinpoint  and  remove  the  different  numerical 
artifacts  such  as  wrap-around  and  ringing. 


8.3.  SAFARI-FIPP  EXAMPLES 

The  use  of  SAFARI-FIPP  for  calculating  the  full  time  response  will  here  be  illustrated 
by  a  few  examples,  all  treating  propagation  in  the  sample  environment  used  for  the 
SAFARI-FIPR  and  SAFARI-FIP  test  cases,  i.e.  a  shallow  water  environment  with  100  m 
water  depth  and  the  layered  bottom  given  in  Table  4. 


■  8.3.1.  SAFARI-FIPP  case  1:  Dispersion  curves 

In  this  and  the  next  example  we  use  SAFARI-FIPP  to  model  the  propagation  of 
seismic  interface  waves  along  the  ocean  sea-bed.  One  of  the  numerical  problems 
particular  to  the  pulse  model  is  the  determination  of  the  necessary  time  window 
in  order  to  avoid  wrap-around.  As  described  above,  the  starting  time  is  easily 
determined  from  the  maximum  wave  speed  in  the  problem.  The  length  of  the  time 
window,  however,  requires  a  little  more  skill  to  determine,  in  particular  in  cases 
where  very  slow  interface  waves  are  present.  However,  SAFARI-FIPP  has  the  option 
of  calculating  group  velocities  of  the  interface  wave,  which  is  a  very  valuable  tool  in 
this  regard.  In  this  first  example  we  will  therefore  calculate  the  dispersion  curves  in 
the  frequency  interval  1-12  Hs  for  the  slowest  fundamental  interface  mode  associated 
with  the  environment  treated  in  SAFARI-FIP  case  1.  To  do  this,  the  data  file  given 
in  Table  18  is  set  up. 

The  resulting  dispersion  curves  are  shown  in  Fig.  20.  The  minimum  group  velocity 
is  observed  to  be  280  m/s.  If  the  pulse  response  has  to  be  calculated  out  to  a  range 
of  2.6  km  and  no  reduced  time  is  applied,  the  length  of  the  time  window  required  in 
order  to  avoid  wrap-around  is  T  =  10  s. 
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Table  18 

SAFARI-FIPP  case  1  data  file 


Data  file 

Description 

Notes 

SAFARI-FIPP  case 

1 

title 

T  G 

options 

1 

3 

centre  frequency  3  Hz 

2 

4 

number  of  layers 

0  C  0  0  0  0  0 

layer  1  (vacuum  halfspace) 

0  1500  00010 

iayer  2  (water  column) 

100  1600  400  0.2 

0.5 

1.8 

0 

layer  3  (silt  sediment) 

120  1800  600  0.1 

0.2 

2.0 

0 

layer  4  (sand  sub-bottom) 

95 

source  depth  95  m 

100  100  1 

receiver  depth  100  in 

300  1E8 

phase  velocity  interval 

3 

250  1  250  0 

wavenumber  sampling  parameters 

4 

612  1.0  12.0  .01 

0.5 

0.5 

0 

frequency  sampling  parameters 

5 

1 

number  of  modes 

6 

0  12  20  2 

frequency  axis  parameters 

200  600  12  100 

velocity  axis  parameters 

1  It  is  indicated  by  option  G  that  the  dispersion  characteristics  of  one  or  more  modes 
should  be  determined.  Option  V  indicate?  that  the  depth-dependent  Green’s  function 
for  the  vertical  particle  velocity  is  used  for  finding  the  peak  corresponding  to  the 
interface  mode. 

*  The  centre  frequency  is  specified  as  3  Hz,  but  this  parameter  has  no  influence  on  the 
dispersion  calculation. 

3  The  wavenumber  interval  is  chosen  so  large  that  the  fundamental  interface  mode  is 
known  to  be  included  at  the  maximum  frequency  of  12  Hz.  Note  that  option  B  was 
not  specified.  The  wavenumber  spectrum  is  therefore  constant  and  equal  to  the  one 
determined  by  the  maximum  frequency. 

*  The  wavenumber  sampling  is  not  extremely  critical  for  the  G  option  because  a  bisection 
technique  is  used  to  iterate  onto  the  correct  position  of  the  modal  peak. 

'  The  frequency  interval  has  been  chosen  as  1  to  12  Hz,  and  to  obtain  a  frequency 
sampling  interval  of  A /  =  0.2  Hz,  the  time  sampling  parameters  have  been  adjusted 
to  yield  a  total  time  window  of  length  T  =  1/A/  ~  5  s.  Since  NRAN  =  0,  no  pulse 
plots  are  generated,  but  we  have  in  any  case  chosen  At  =  0.01  a,  which  is  the  sampling 
required  to  nicely  represent  graphically  a  wave  of  frequency  12  Hz.  The  appropriate 
number  of  time  samples  therefore  follows  as  NT  =  29  =  512. 

6  In  order  to  determin  e  the  width  of  the  time  window  for  the  following  pulse  calculations, 
only  the  group  velocity  of  the  fundamental  interface  mode  is  required. 
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DISPERSION  CURVES 


Frequency  (Hz) 


Fig.  20.  SAFARI-FIPP  case  1:  dispersion  curves  for  the  fundamental  in¬ 
terface  mode  in  a  stratified  sea-bed.  The  solid  curve  indicates  the  phase 
velocity  and  the  dashed  curve  the  group  velocity. 
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■  8.3.2.  SAFARI-FIPP  case  2:  Seismic  interface  wave  propagation 

In  this  example  we  calculate  the  time  response  of  5  vertical  geophones  on  the  ocean 
bottom,  spaced  500  m  apart  out  to  a  range  of  2.5  km,  using  a  source  of  type  2  with 
a  centre  frequency  of  5  Hz.  The  environment  is  exactly  as  in  the  previous  example. 
This  is  done  by  specifying  the  data  file  given  in  Table  19. 

The  resulting  stacked  plots  of  the  time  behaviour  of  the  vertical  ai.d  horizontal  par¬ 
ticle  velocities  are  shown  in  Fig.  21a  and  b,  respectively.  Note  that  the  amplitudes 
have  been  multiplied  by  range  before  being  plotted.  The  highly  dispersive  nature  of 
the  slow  seismic  interface  wave  is  characteristic  for  this  type  of  propagation  situa¬ 
tion.  Although  not  particularly  evident  in  the  present  case,  it  is  also  characteristic 
that  the  horizontal  component  shows  a  second  interface  mode  more  pronounced  than 
the  vertical  component.  For  a  more  detailed  description  of  the  behaviour  of  seismic 
interface  waves,  reference  is  made  to  [16-18].  In  Fig.  21b  for  the  horizontal  compo¬ 
nent,  the  first  arrival  is  the  lowest  waterborne  normal  mode,  which  is  only  weakly 
present  on  the  vertical  component.  This  arrival  also  contains  a  head  wave  arising 
from  the  sub-bottom,  but  this  wave  is  not  easily  isolated  from  the  waterborne  mode. 
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Table  19 

SAFARI-FIPP  case  2  data  file 


Data  file 

Description 

Notes 

SiFARI'FIPP  case  2 

title 

f  H  S  2  F  J 

options 

1 

5  0 

centre  frequency  5  Hz 

2 

4 

number  of  layers 

0  0  0  0  0  0  0 

layer  1  (vacuum  halfspace) 

0  1600  00010 

layer  2  (water  column) 

100  1600  400  0.2  0.S  1.8  0 

iayer  3  (silt  sediment) 

120  1800  600  0.1  0.2  2.0  0 

layer  4  (sand  sub-bottom) 

96 

source  depth  95  m 

100  100  1 

receiver  depth  100  m 

300  1E8 

phase  velocity  interval 

3 

1024  1  960  0 

wavenumber  sampling  parameters 

4 

2048  0.0  12.6  .006  0.6  0.6  5 

frequency  sampling  parameters 

5 

0 

reduction  velocity 

6 

0  10  20  2 

time  axis  for  stacked  plots 

7 

0  3  16  0.6 

range  axis  for  stacked  plots 

8 

1  Options  V  and  E  are  specified  in  order  to  cal  ulate  both  vertical  and  horizontal  particle 
velocities.  Since  option  S  has  been  selected,  the  time  responses  will  be  displayed  in  a 
range-stacked  format.  A  source  pulse  of  type  2  has  been  selected.  Option  F  invokes 
the  Filon  integration  scheme  which  requires  less  sampling  points  than  the  default 
trapesoidal-rule  integration  scheme.  Finally,  the  J  option  is  specified  to  move  the 
wavenumber  integration  contour  out  into  the  complex  plane.  This  is  not  necessary 
for  the  interface  wave,  which  is  itself  highly  attenuated,  but,  as  was  demonstrated  in 
SAFARI-FIP  case  1,  normal  modes  will  become  important  for  frequencies  above  the 
centre  frequency  of  5  Hs.  Therefore  a  significantly  larger  number  of  sampling  points 
than  the  actual  1024  would  be  required  at  the  higher  frequencies  in  the  band.  Note 
that  option  B  was  not  specified  because  the  exponential  decay  of  the  wavenumbi. 
integration  kernel  will  be  very  slow  at  the  low  frequencies  in  cases  where  interface 
waves  are  important. 

1  The  centre  frequency  is  specified  as  5  Hz.  By  inspection  of  Fig.  16b  it  is  clear  that 
the  source  pulse  will  contain  significant  energy  up  to  a  frequency  of  12. &  Hs,  which  is 
therefore  selected  as  the  upper  limit  for  the  calculations.  The  default  contour  offset  is 
applied  by  specifying  COFF  =  0. 

3  As  in  case  1,  a  wide  wavenumber  spectrum  is  used  in  order  to  encompass  all  real  angles 
of  propagation  and  all  evanescent  waves  that  are  present. 

4  The  number  of  wavenumber  samples  selected  is  1024,  which  translates  into  Afc  rBtl  ~ 
0.64  <  i.e.  sufficient  also  for  the  standard  trapezoidal  rule  integration  and  therefore 
more  than  enough  for  the  selected  Filon  scheme.  The  last  parameter  IIMC  is  set  to  0 
in  order  not  to  create  any  wavenumber  integrand  plots.  By  selecting  XC2  as  950  the 
tapering  is  invoked,  to  take  care  of  possible  discontinuities  at  the  high  wavenumber 
truncation  point. 
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The  frequency  interval  selected  is  0  to  12.5  Hi,  which  L-  an  appropriate  choice  for 
pulse  type  2,  see  Fig.  21.  The  total  time  window  is  chosen  a  little  larger  than  10  s  by 
specifying  DT  =  0.006  and  IT  =  2048.  The  ranges  where  pulses  are  to  be  calculated  are 
determined  by  the  iast  3  parameters.  In  the  present  case  there  are  5  receiver  ranges 
spaced  500  m  apart  in  the  interval  500  to  2500  m.  Note  that  ranges  are  always  given 
in  km. 

Due  to  the  short  ranges  involved,  no  time  reduction  is  applied;  therefore  CRID  =  0. 
The  time  axis  on  all  generated  plots  is  specified  to  contain  a  window  of  width  10  s, 
starting  at  0  s.  The  length  of  the  time  axis  selected  is  20  cm. 

Since  the  receivers  are  placed  at  ranges  between  0.5  and  2.5  km,  the  range  axis  for  the 
stacked  plots  has  been  chosen  to  range  from  0  to  3.0  km  in  order  not  to  truncate  the 
first  and  last  trace. 
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■  8.3.3.  SAFARI-FIPP  case  3:  Normal  mode  propagation 

In  this  last  example  the  stacking  in  depth  instead  of  range  will  be  demonstrated  by 
calculating  the  field  produced  at  a  vertical  array  of  hydrophones  for  a  source  at  50  m 
depth  and  with  range  offsets  of  5  and  10  km.  The  environment  is  the  same  as  treated 
in  examples  SAFARI-FIP  2  and  3,  and  the  source  is  assumed  to  be  of  type  4  with  a 
centre  frequency  of  30  Hs.  Only  the  discrete  part  of  the  wavenumber  spectrum  wil 
be  included  in  order  to  demonstrate  how  the  truncation  of  the  integration  interval 
is  introduced  in  the  pulse  calculations.  The  data  file  is  given  in  Table  20. 

The  resulting  stacked  plots  are  shown  in  Fig.  22.  Note  the  two  arrivals  corresponding 
to  the  two  normal  modes.  As  expected,  the  first  arrival  is  the  first  mode  with  its 
largest  amplitude  at  approximately  60  m  depth.  The  second  mode  arrives  0.2  to  0.3  s 
later,  and  it  has  a  low  anplitude  at  60  rn  depth.  By  comparing  the  results  for  the 
two  different  ranges,  the  characteristic  dispersion  of  the  modes  is  evident. 
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SAFARI-FIPP  c»«i  3  data  file 


Data  file 

Description 

Notes 

SAFttl-FIPP  case  3 

title 

N  0  4  J  F  B 

options 

1 

30  0 

centre  frequency  30  Hs 

6 

number  of  layers 

0  0  0  0  0  0  0 

layer  1  (vacuum  halfspace) 

0  1600  -1480  0  0  t  0 

layer  2  (water  layer) 

30  1480  -1490  0010 

layer  3  (water  layer) 

100  1800  400  0.2  0.5  1.8  0 

layer  4  (silt  sediment) 

120  1800  600  0.1  0.2  2.0  0 

layer  5  (sand  sub-bottom) 

60 

source  depth  50  m 

90  10  9 

receiver  depths  10  to  90  m 

2 

1460  1700 

phase  velocities  CHIN  and  CNAX 

3 

650  60  600  0 

wavenumber  sampling  parameters 

3 

1024  16  75  0.0012  552 

frequency  sampling  parameters 

4 

1600 

reduction  velocity 

5 

0  1  20  0.2 

time  axis  for  stacked  plots 

6 

100  0  15  20 

depth  axis  for  stacked  plots 

7 

1  Option  M  indicates  that  the  normal  stress,  equal  to  the  negative  value  of  the  acoustic 
pressure,  should  be  calculated.  The  depth  stacking  is  invoked  by  specifying  option 
D,  and  the  wavenumber  contour  offset,  option  3,  is  again  applied.  Option  F  selects 
the  Filon  integration  scheme,  and  finally  option  B  indicates  that  the  phase  velocity 
interval  (or  slowness  interval)  should  be  held  constant  at  all  frequencies  in  contrast  to 
the  former  example,  where  the  wavenumber  interval  was  he.d  constant.  The  B  option 
is  essential  when  the  calculations  involve  only  a  limited  beamwidth. 

3  There  are  9  receivers  in  the  vertical  array,  with  the  lowermost  receiver  at  90  m  d  ii 
and  the  uppermost  at  10  m  depth.  The  deepest  receiver  has  been  specified  first,  since 
the  amplitude  scaling  is  controlled  by  the  first  receiver  depth,  and  the  amplitude  is 
expected  to  be  larger  at  the  bottom  than  at  the  pressure-release  surface. 

3  As  in  e^arrr'  "AFARI-FIP  case  3,  the  wavenumber  interval  has  been  selected  to 
include  011’*'  •  j  propagating  normal  modes,  and  the  tapering  is  applied  by  properly 
specifying  XC1  and  IC2.  No  integrand  plots  have  been  requested  (IIMC  as  0),  although 
this  feature  is  often  applied  to  check  the  wavenumber  truncation. 

*  As  can  be  observed  from  Fig.  18b,  the  selected  frequency  interval  15  to  75  Hs  will 
include  the  whole  main  lobe  of  the  source  frequency  spectrum.  The  time  sampling 
parameters  have  been  chosen  to  yield  a  total  time  window  of  1.2  s,  which  turns  out  to 
be  sufficient  for  the  ranges  selected.  As  above,  a  dispersion  calculation  for  the  normal 
modes  (2  at  the  centre  frequency)  could  have  been  performed  in  advance. 

1  Since  the  fast-bottom  head  wave  hes  been  excluded  from  the  wavenumber  integration 
interval,  no  arrivals  fa**'  \u»  the  water  velocity  (1500  m/s)  are  expected.  Therefore 

this  velo'-:‘”  '  <ten  '  as  reduction  velocity. 

8  The  time  uua  covers  .  .e  first  1  s  of  the  selected  time  window. 

T  The  depth  axis  for  the  stacked  plots  covers  the  whole  water  column,  leaving  enough 
space  for  both  the  first  and  the  last  trace. 
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NORMAL  STRESS 


NORMAL  STRESS 


Fig.  22.  Results  for  SAFARI-FIPP  case  3:  depth-stacked  plots  of  hy¬ 
drophone  traces  at  two  different  ranges:  (a)  5  km;  (b)  10  km. 
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Appendix  A 

Installing  and  running  SAFARI  on  the  FPS164 

SAFARI  is  installed  at  SACLANTCEN  on  both  the  VAX  8600  cluster  and  the  FPS164 
Attached  Processor.  The  installation  on  the  VAX  proceeds  exactly  as  described  in 
Sect.  5.  For  the  FPS164  the  simulation  of  the  FPS  library  functions  should  of  course 
not  be  used.  Further,  the  most  critical  subroutines  have  been  written  in  APAL64 
assembler  and  a  special  set  of  subroutines  is  available  for  doing  asynchronous  I/O 
to  all  scratch  files.  This  last  feature  requires  the  FPS164  to  be  equipped  with  a  D64 
disk  subsystem  and  to  run  under  the  Single  Job  Executive  (SJE)  operating  system. 

The  FPS1S4  version  is  compiled  and  linked  by  means  of  the  following  command  file: 


COMFPS.COM _ 

$  SET  DEF  DS5 : [SCHNI . ANIS] 

$  ! 

$  !  SUBROUTINES 
i  i 

$  APFTN84/X0FF-ALL/0PT-2  FIPASJE30 
I  APFTN64/I0FF«ALL/0PT»2  FIPSSJE30 
$  APFTN64/X0FF-ALL/0PT-2  FIPPSJE30 
3  APFTNB4/X0FF-ALL/0PT-2  FIPUSJE30 
$  APFTN64/X0FF»ALL/0PT«2  FIPTSJE30 
$  APFTN«4/X0FF»ALL/0PT*2  FIPRSJE30 
$  ! 

t  !  ASSEMBLER  SUBROUTINES 

I  I 

$  APALS4  APCBGE5 
I  APALS4  CYMOVI 
$  APAL64  CYIM0Y 
$  ! 

$  !  SUBROUTINES  FOR  ASYNCHRONOUS  I/O 

$  ! 

)  APFTN64/X0FF»ALL/0PT«2  ASI04 
$  ! 

$  !  CREATE  LIBRARY 

$  ! 

$  APLIBR«4/0UT»SAFLIB/IN«(FIPPSJE30,FIPSSJE30,FIPUSJB30,- 

FIPTS JE30 , FIPRS JE30 .FIPAS JE30 , APCBGE5 , CYMOVI , CTIHOt , ASI04) 

$  ! 

I  I  MAIN  PROGRAMS 

I  I 
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I  APFTM64/X0FF-ALL  SJEFIPP30 
$  APLINK64  S JEFIPP30 , SAFLIB . AOL 
I  APFTN64/X0FF*ALL  SJEFIP30 
I  APLINK64  S JEFIP30 , SAFLIB . AOL 
I  APFTN64/X0FF=ALL  SJEFIPR30 
I  APLINK64  SJEFIPR30, SAFLIB. AOL 
$  ! 

9  !  CREATE  AND  SUBMIT  COMMAND  FILE  FOR  COPTINC 
9  !  IMAGES  INTO  FPS164/SJE  FILES,  HERE  UNDER 
9  !  ACCOUNT  :SCHMI 

$  ! 

9  OPEN/WRITE  COM  COP.COM 
9  WRITE  COM  "I  SET  DEF  US5 : LSCHMI . ANIS] " 

$  WRITE  COM  "$  SJE" 

$  WRITE  COM  "ATT/W" 

9  WRITE  COM  "ACC  :SrHMI" 

9  WRITE  COM  "COPTIN/BI  SJEFIP30.IKG.SJEFIP30" 

9  WRITE  COM  "COPTIN/BI  SJEFIPP30.IKG,SJEFIFP30" 
9  WRITE  COM  "COPTIN/BI  S JEFIPR30 . IMG , S JEFIPR30" 
9  WRITE  COM  "DEI" 

9  WRITE  COM  "QUIT" 

9  CLOSE  COM 
9  FPSqUE  COP 


Also,  a  special  command  file  is  required  for  running  SAFARI  on  the  FPS164,  e.g.  for 
the  transmission  loss  program  SAFARI-FIP: 


FIP.COM 


SET  DEF  US6: [SCHMI.ANIS] 

0PEN/W1 ITE  COM  SAFARIFIP.COM 
WRITE  COM  "9  SJE" 

WRITE  COM  "ATTACH/W  1" 

WRITE  COM  "COPTIN  ’  ’PI ’ .DAT.FTNOOl" 
WRITE  COM  " : SCHMI : S JEFIP30" 

WRITE  COM  "COPTOUT  FTN019,  "PI ’ .PLP" 
WRITE  COM  "COPTOUT  FTN020,  "PI » .PLT" 
WRITE  COM  "COPTOUT  FTN028,  "PI ’ .CDR” 
WRITE  COM  "COPTOUT  FTN029,  "PI » .BDR” 
WRITE  COM  "DETACH" 

WRITE  COM  "QUIT" 

CLOSE  CON 
•SAFARIFIP 

ASSIGN/USER  ’PI’. PLP  FOR019 
ASSIGN/USER  ’PI ' .PLT  F0R020 
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$  ASSIGN/9SER  SYSSCONMAND :  SYSSINPUT 
9  RUN  USE:  |[SCHMI . YAXFIP] FIPPLOT 
$  ASS/USER  ’PJ’.CDR  F0R055 
$  AS3/USER  ’Pl'.BDR  F0R017 
!)  ASS/USER  STSICOMMAND :  SYSSINPUT 
$  RUN  USE:  [SCHMI.FIPC0N1C0NTUR.EXE 

The  two  other  modules  are  run  by  means  of  similar  command  files. 
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Appendix  B 
Running  FIPPLOT 

When  SAFARI  is  run  by  using  the  command  files  given  in  Sect.  5  and  Appendix  A, 
the  plot  program  FIPPLOT  is  automatically  executed  to  obtain  the  graphical  output. 
The  plot  data  are  transferred  through  two  files.  The  file  with  extension  PLP  contains 
plot  parameters  such  as  axis  lengths,  titles,  labels  etc,  whereas  the  other  Ale,  with 
extension  PLT,  contains  the  actual  data  values  to  be  plotted.  Both  files  are  formatted 
ASCII  files,  and  it  is  therefore  possible  to  change  the  layout  of  the  plots  by  editing 
the  PLP  file  and  then  execute  FIPPLOT  independently.  Further,  FIPPLOT  has  such 
a  general  setup  that  it  is  easily  interfaced  to  any  other  numerical  code  requiring 
graphical  output.  We  will  therefore  describe  here  the  structure  of  the  PLP  file  in 
detail. 

As  a  characteristic  example,  the  PLP  file  generating  the  plot  shown  in  Fig.  11a  is  as 
follows: 


u1024 

uFIPPuuSTLDiV , CPX , IT A 

UDEPTH  AVERAGED  LOSS 
uSAFARI-FIP  case  3. 

modulo 

u2 

uFr«q:  30.0  Hz$ 

USD :  SO . 0  at 

number  of  labels 

u20. 000000 

ELEN 

u12. 000000 

TLEX 

uO 

grid  type 

u0. 000000 

XLEFT 

u6. 00000 

X11GHT 

ui. ooooo 

xmc 

ui. ooooo 

uRanga  (ka)$ 

ULIN 

XDIV 

u80.0000 

TD0VN 

u20.0000 

TOP 

ulO.OOOO 

TI1IC 

u 1.00000 

uXoraal  stress  (dB//lPa)t 
uLIM 

TDIT 

ul 

HC 

yll2 

N 

u0. 044967 

XHIX 

u0. 449667 -001 

DX 
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THIN 

nr 


u0. 000000 
u0. 000000 
yFIPyyyPLlKND 


Notice  especially  the  initial  space  (denoted  by  u)  which  is  necessary  for  FIPPLOT  to 
correctly  read  the  file.  The  PIP  file  will  always  start  with  a  block  size  modulo,  with 
which  the  corresponding  data  hie  was  written.  This  parameter  is  a  leftover  from 
the  versions  using  binary  data  files,  but  for  the  present  version  of  FIPPLOT  this 
parameter  is  dummy.  The  file  ends  with  the  PLTEND  flag  signalling  the  end-of-file. 

Between  these  two  records  the  blocks  of  actual  plot  parameters  are  specified,  in 
the  present  case  for  a  single  plot  only.  Any  number  of  blocks  could,  however,  be 
included,  each  generating  one  plot. 

The  plot  parameter  block  starts  with  a  record  specifying  a  12-character  plot  iden¬ 
tification  (FIPPuySTLDAV)  followed  by  a  series  of  3-dharactcr  options  separated  by 
commas.  The  files  generated  by  the  SAFARI  modules  will  in  general  not  have  any 
options  specified,  but  these  are  important  tools  for  generating  final  plots.  The  fol¬ 
lowing  options  are  currently  available: 

DUP:  The  duplex  character  generator  will  be  used  in  stead  of  the  default  simplex. 

CPX:  The  complex  character  generator  is  selected.  This  is  the  option  chosen  for 
all  plots  in  this  report . 

ITA:  The  italic  character  generator  is  used. 

IXA:  Integer  format  will  be  uned  for  plotting  the  z-axis  tick  mark  numbers  instead 
of  the  default  decimal  format. 

ITA:  Integer  format  will  be  used  for  plci  ting  the  y-axis  tick  mark  numbers  instead 
of  the  default  decimal  format. 

DSD:  If  the  plot  contains  more  than  one  curve  (SIC  >  1),  then  the  first  curve  will 
be  plotted  with  a  solid  line,  the  second  with  a  dashed  line  and  the  third 
with  a  dotted  line.  If  more  than  three  curves  are  plotted,  this  sequence  will 
be  repeated.  To  generate  the  dispersion  plot  in  Fig.  20  the  DSD  option  was 
applied. 

MRK:  A  marker  will  be  plotted  for  ea^h.  10th  data  point  on  the  curves.  In  the  case 
of  more  curves,  different  markers  will  be  used  for  each. 

COL:  This  option  will  plot  individual  curves  in  different  colours,  using  the  repeat- 
able  sequence:  red-green-blue-cy an-magent a-yellow . 

SCA:  The  default  character  size  has  been  selected  to  yield  a  reasonable  ‘balance’ 
of  the  plot  when  the  length  of  the  z-axis  is  20  cm.  This  option  will  scale 
the  character  size  accordingly  for  larger  or  smaller  sizes.  This  is  particularly 
important  when  preparing  illustrations  for  reports  and  journal  papers,  where 
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a  uniform  character  size  is  desirable,  independent  of  the  size  of  the  original 
plots. 

RES:  When  this  option  is  selected,  FIPPLOT  will  request  a  scaling  factor  (1.0  = 
no  scaling)  which  will  be  applied  to  all  curves  in  the  plot,  but  not  to  the  axis 
values.  This  option  is  particularly  useful  for  the  range  and  depth-stacked 
tune  series  produced  by  SAFARI-FIPP. 

TCT:  This  option  is  used  in  connection  with  the  stacked  time  series  plots  in  order 
to  truncate  the  amplitude  of  each  trace  at  a  value  corresponding  to  half  the 
distance  between  traces.  Used  in  particular  to  avoid  overlap  of  curves  when 
rescaling  with  a  large  factor  (option  RES). 

NOP:  With  this  option  both  the  parameters  from  the  PLP  file  and  the  data  from 
the  PLT  file  are  read  but  no  plot  is  produced.  It  is  therefore  used  for  time¬ 
saving  when  only  selected  plots  are  required. 

The  next  record  specifies  the  main  title  of  the  plot,  followed  by  a  record  containing 
the  title  specified  in  the  SAFARI  data  file  for  the  actual  run.  This  title  will  be  plotted 
just  above  the  plot  frame. 

Then  follows  a  sub-block  containing  the  labels  to  be  plotted  in  the  upper  right  corner 
of  the  plot  frame.  The  number  of  labels  (>  0)  is  given  first,  followed  by  the  label 
texts,  each  of  which  should  be  on  separate  lines  and  terminated  by  a  $. 

The  parameters  XLEN  and  YLEN  specify  the  length  in  cm  of  the  z-  and  y-axis,  re¬ 
spectively.  The  parameter  labelled  ‘grid  type’  indicates  whether  a  grid  should  be 
plotted.  A  value  of  1  will  produce  a  dotted  grid  as  in  Pig.  20. 

The  next  6  records  contain  the  parameters  for  the  z-axis  of  the  plot.  XLEFT  and 
XRItiHT  are  the  data  values  at  the  left  and  right  borders  of  the  plot  frame,  respec¬ 
tively,  whereas  XINC  is  the  distance  in  the  same  units  between  the  tick  marks.  XDIV 
is  a  multiplication  factor  which  will  be  applied  to  both  the  axis  parameters  and  the 
data  values.  After  XDIV  the  z-axis  label  is  specified,  terminated  by  a  $,  and  finally 
LIN  indicates  that  the  z-axis  should  be  linear.  Another  possibility  is  a  logarithmic 
axis,  which  has  not  yet  been  implemented,  however.  The  parameters  for  the  y-axis 
are  given  in  the  same  way  in  the  next  6  records. 

The  parameter  NC  specifies  the  number  of  curves  to  be  plotted.  For  each  curve 
a  sub-block  of  5  records  has  to  be  specified.  The  first  parameter  N  indicates  the 
number  of  points  in  the  curve.  If  N  is  negative,  no  curve  will  be  plotted;  instead  a 
marker  will  be  plotted  at  the  position  of  each  data  point.  The  parameter  XMIN  is  the 
z-coordinate  (range)  of  the  first  data  point,  whereas  DX  is  the  equidistant  spacing. 
If  DX  had  been  specified  as  0,  then  the  N  z-coordinates  of  the  data  points  would  be 
read  from  the  PLT  file.  In  that  case  XMIN  would  be  interpreted  as  an  z-offset  to  be 
applied  to  the  curve.  The  same  rules  apply  to  THIN  and  DY.  In  the  actual  case  the 
y- values  (transmission  losses)  will  therefore  be  read  from  the  PLT  file,  and  no  y-offset 
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will  be  applied.  The  offsets  are  only  used  for  the  stacked  time  series  plots,  where 
it  is  important  to  be  able  to  rescale  (see  option  RES  above)  the  amplitudes  without 
changing  the  trace  offset.  If  both  DZ  and  DT  are  specified  as  0,  then  FIPPLOT  will 
first  read  all  N  i- values  and  then  all  y- values.  As  an  example  of  this,  reference  is 
made  to  the  PLP  file  for  generating  the  sound-speed  profile  plot  inserted  in  Fig.  10a. 

When  the  PLP  file  has  been  edited,  FIPPLOT  is  executed  with  the  following  com¬ 
mand  file: 


$  ASS/USER  ’PI* .PLP  F0RO19 
$  ASS/USER  ’Pl’.PLT  FOR020 
$  ASS/USER  STSICOMMAND :  STSRINPUT 
$  RUN  FIl’PLOT 


All  plots  presented  in  this  report  have  been  edited,  and  they  therefore  do  not  appear 
exactly  as  those  produced  automatically. 
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Appendix  C 
Running  CONTUR 

As  was  the  case  for  FIPPLOT,  the  command  files  given  in  Sect.  5  and  Appendix  A  will 
automatically  execute  the  contour  plotting  program  CONTUR.  The  plot  parameters 
are  transferred  from  SAFARI  into  the  file  with  extension  CDR,  whereas  the  actual 
contour  data  are  placed  in  the  file  with  extension  BDR.  The  data  file  is  allowed  to 
be  in  both  ASCII  or  binary  format,  but  as  SAFARI  is  run  both  on  the  VAX  and  the 
FPS164  at  SACLANTCEN,  it  has  been  most  convenient  to  use  the  ASCII  format. 

The  CDR  file  can  be  edited  to  change  the  layout  of  the  contour  plot.  As  a  typical 
example,  the  CDR  file  used  to  create  the  contour  plot  of  Fig.  14b  is  shown  below: 


CONOR , FIP , FHT , CPX , ONI , TEK 
SAFARI -FIP  case  6. 

SAFFIP6 .BDR; 1 
Rang*  (■) 

0.0000 
299.8535 
0.0000 
300.0000 
16.0000 
60.0000 
Depth  (a) 

60.0000 

125.0000 

6.2500 

25.0000 

141.0000 

61.0000 

1.0000 

1.0000 

0.0000 

60.0000 

126.0000 

60.0000 

141.0000 

61.0000 

1000.0000 

0.0000 

6.0000 

6.0000 

21.0000 


RHXN 

RJUI 

XLEFT 

XRIGHT 

XSCALE 

XINC 

TOP 

TD0VN 

TSCALE 

TINC 

data  points  along  z-axis 

data  points  along  y-axis 

DITX 

DIVT 

FLAGRC 

RDUP 

RDL0 

source  depth  (m) 

grid  points  along  z-axis 

grid  points  along  y-axis 

frequency  (Hs) 

dummy 

CAT 

NANG 

ZHIN 
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54.0000 

Z  MX 

3.0000 

ZINC 

2.0000 

z-origin  of  plot  in  inches 

0.0000 

dummy 

2.0000 

y-origin  of  plot  in  inches 

0.0000 

NSN 

0.1000 

HCTPT 

0.1400 

HGTC 

-3.3000 

LABPT 

1.0000 

NDIV 

5.0000 

NARC 

-1.0000 

LABC 

-1.0000 

LVGT 

The  first  record  specifies  one  5-character  option  (CONOR)  followed  by  a  series  of  3- 
letter  options.  The  CONDR  option  indicates  that  the  actual  contour  plot  is  of  the 
depth-range  type  and  CONTUR  will  interpret  the  parameters  accordingly.  This 
option  should  therefore  never  be  changed.  The  first  3-letter  option  (FIP)  is  purely 
for  identification  and  has  no  further  effect.  The  FMT  option  indicates  that  the  BDR 
data  file  is  ASCII  formatted  (BIN  for  binary  format).  These  first  3  options  should 
always  be  present  in  the  specified  order,  but  the  options  following  are  optional  and 
can  be  given  in  any  order.  The  implemented  options  are  as  follows: 

DUP:  The  duplex  character  generator  will  be  used  for  DISSPLA  plots  (MINDIS  at 
SACLANTCEN )  instead  of  the  default  simplex. 

CPX:  The  complex  character  generator  is  selected. 

ITA:  The  italic  character  generator  will  be  used. 

UNI:  The  UNIRAS  plot  package  will  be  used  to  generate  a  colour  or  grey-tone 
contour  plot.  DISSPLA  (MINDIS  at  SACLANTCEN)  is  the  default  plot 
package. 

VTT:  The  raster  plot  will  be  produced  on  a  VT240  terminal.  This  only  has  effect 
if  option  UNI  is  specified. 

G41:  The  raster  plot  will  be  produced  on  a  Tektronix  4105  terminal.  This  only 
has  effect  if  option  UNI  is  specified. 

TEK:  The  raster  plot  will  be  produced  on  a  Tektronix  4691  ink  jet  plotter.  This 
only  has  effect  if  option  UNI  is  specified. 

PRX:  The  raster  plot  will  be  produced  on  a  Printronix  printer.  This  only  has 
effect  if  option  UNI  is  specified. 

COL:  If  the  option  UNI  has  been  selected,  a  colour  raster  plot  using  au  in-built 
colour  scale  will  be  generated.  By  default  a  grey-tone  scale  is  applied. 
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VUG:  If  option  UNI  is  selected,  a  viewgraph  will  be  produced  on  a  Tektronix  4691 
or  4692  ink  jet  plotter. 

ROT:  The  UNIRAS  plot  will  be  rotated  90°.  Mainly  used  for  generating  viewgraphs 
of  plots  with  the  z-axis  being  larger  than  the  y-axis. 

After  the  option  record  there  is  a  record  containing  the  title  of  the  plot  and  a  record 
containing  the  name  of  the  file  containing  the  data,  i.e.  the  BDR  file.  Except  for 
2  records  containing  the  z-axis  and  y-axis  labels,  the  rest  of  the  records  contain 
numerical  parameters,  all  supplied  with  a  descriptive  label.  In  general  only  i  few  of 
these  parameters  should  be  changed.  The  most  important  ones  are  described  below. 

The  lengths  of  the  z-  and  y-axes  are  controlled  by  the  parameters  XSCALE  and 
YSCALE,  respectively.  CONTUR  requires  these  parameters  to  be  specified  in  coor¬ 
dinate  units  per  cm.  ZHIN  and  ZNAX  specify  the  limits  of  the  contouring  interval, 
whereas  ZINC  is  the  contour  increment.  If  UNIRAS  is  selected,  areas  with  small  data 
values  will  be  coloured  with  magenta  (black  in  the  grey-tone  mode).  The  colour  scale 
then  moves  through  different  red  tones  into  blue  (white  in  the  grey-tone  mode).  If, 
however,  ZMIN  >  ZNAX,  then  the  colour  scale  will  be  inverted  by  CONTUR.  If  the 
default  DISSPLA  package  is  selected,  only  contour  lines  with  identifying  numbers 
are  plotted.  In  the  present  case,  removal  of  the  UNI  option  leads  to  the  contour  plot 
shown  in  Fig.  14a. 

The  last  important  parameter  is  NSN,  which  controls  the  amount  of  smoothing  ap¬ 
plied  to  the  calculated  contours.  This  parameter  can  be  set  to  any  value  between 
0  and  10,  with  0  corresponding  to  no  smoothing.  It  is  obvious  that  this  parameter 
should  be  used  with  extreme  care. 

After  editing  the  CDR  file,  CONTUR  may  be  executed  by  means  of  the  following 
command  file: 


$  ASS/DSER  ’PI ’ . CDR  F0R065 
$  ASS/DSER  ’PI ’.BDR  F0R017 
$  ASS/OSKR  SYStCOHMAND:  STSIINPOT 
$  RUN  OSS : [SCHHI . FIPCOi] CONT0R 
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